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DA&OS  MocUfr-Lcation 

Initiator:  Mrs.  I.  Hussey 

Problem  No:  3018-3  Project  No:  0001 

The  purpose  of  this  work  was  to  combine  into  one  program  two  earlier 
versions  of  DABOS,  which  is  an  orbit  determination  and  ephemeris  computation 
program.  From  a mathematical  point  of  view,  the  difference  between  the  two 
earlier  versions  was  that  one  used  the  Jacchia  1964  atmospheric  density  model 
and  the  other  used  the  Standard  Atmosphere  Supplements  1966.  These  two  density 
models  are  rather  similar  and  it  was  possible  to  intertwine  the  formulae  in  a 
very  efficient  manner,  allowing  the  user  to  select  which  model  was  to  be  used 
for  a particular  run. 

Shotit-Tvm  PoiAieji  Spectra  otf  TAopo4ca£teA  Slgnati 

Initiator:  Mr.  U.  Lammers 

Problem  No:  30  3 Project  No:  8682 

This  problem  required  the  analysis  of  sets  of  equally  spaced  data  (X(I)). 

A program  was  written  which  transformed  the  original  input  data  (X(I))  as 
follows: 

Y(I)  = lO*^20 

The  resulting  Yflj's  were  then  Fourier  analyzed  using  a Fast  Fourier 
Transform  and  the  results  graphically  displayed  for  the  Initiator.  In 
addition,  the  program  was  written  with  the  option  of  analyzing  various  seg- 
ments of  the  data  and  producing  graphic  displays  of  each  segment. 


An algAii  o{,  Electromagnetic  Wave*  and  Turbulence 

Initiator:  Dr.  R.  Fante 

Problem  No:  3026-7  Project  No:  5635 

This  problem  required  the  generation  of  four  computer  programs,  each  of 
which  studied  the  behavior  of  electromagnetic  waves  in  the  ionosphere.  The 
documentation  titles  for  these  four  programs  are: 

1.  Calculation  of  the  Frequency  Spectrum  of  a Laser  Beam  Propagating 
in  Turbulence  with  a Von  (Carman' s Spectrum. 

2.  The  Frequency  Spectrum  of  a Plane  Wave  Propagating  in  Turbulence 
with  Constant  Winds. 

3.  The  Frequency  Spectrum  of  a Plane  Wave  Propagating  in  Turbulence 
with  Random  Winds,  and 

4.  The  Frequency  Spectrum  of  a Plane  Wave  Propagating  in  Turbulence 
with  Random  Winds  varying  along  the  Propagation  Path. 

The  first  program  approximated  SQ(ft)  where: 

SoCflj  dRo  cos(flRo)  [FUfRojJ 

where : 

FU(Ro)  = 2S  * (SIG*XK0*BL0)2/(4ir2)  £ 9t  ^°°  dX<KX,Ro,9) 

with: 


<J>(X,Ro,0)  - X exp | PX2  ♦ iKRo.X.6)} 
<KRo,X,9)  = / dy  F(y) 


2 


, if  1.5  < Z 


F(y)  = [1.5Z2  - 1 . 86Z5'  3 - .254Zll/3]  if  Z < .7 

= -1  ♦ '489'  , if  .7  < Z < 1.5 

(Z)-3  “ “ 


= -1  + 


.994  Zl/s  e"z 


♦ 1 --Z-4 

< 92  81Z2  2187Z3  ( 


where 


Z = [Ro2  + 2Xy  Ro  cos©  + X2y2] 


SIG  = .78  * XKO2  * CN2  * exp(|  An (BLO) ) 


p = - (S*SIG2/(4tt))  + ( 2 + — § ■)  * T 

2 4XK02*S  (4ttF2) 


+ SIG  * S * T/(2ttF))  (XKO* BLO) 2 


The  integral  \j)(Ro,X,0)  was  evaluated  by  finding  all  the  discontinuities  in 
F(y)  and  using  a Simpson's  rule  type,  numerical  integration  technique,  for 
each  of  the  continuous  line  segments  in  F(y)  for  o < y < T.  The  integrals 
over  0 and  X used  the  technique  described  in  reference  1 . Here  the  X 
integration  was  performed  from  <j>  to  XMAX,  where  XMAX  is  an  input  parameter, 
just  as  S,  XKO,  BLO  and  CN2. 

The  integration  over  Ro  was  performed  from  o to  NARo,  where  ARo  and  N 
are  input  parameters,  where 


n - 2ttj  NMAX 

~ WAX* ARo  for  J-0*1’  •••»  — 2~ 


A New  Adaptive  Simpson  Integration  Routine,  Neil  Grossbard,  Space  Data 
Analysis  Laboratory,  Boston  College,  Chestnut  Hill,  Mass.  02167;  AFCRL-70-0504, 
Scientific  Report  No.  1,  September  1970. 
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Jeo  S0(i<jJ's  wore  calculated  by  using  a NMAX  long  Fast  Fourier  Transform  to 
simulate  a trapezoidal  rule.  Here  NMAX  = 2^  (U,  integer)  and  2*"1  < 2N  < NMAX. 

The  second  program  approximated  S (ft)  where: 

S ^ (ft)  = ^ fQ  du  cos(ftu) 


where 


i + 1.5  n2  - 1.86  ns/3  - 0.254  n 3 if  0 < n < .7 


.489 


if  .7  < n < 1.5 


7 . 175  l .c  ^ , r 

+ > if  n > i.5 

8in2  2i87n3  ) 


As  in  the  first  program  a set  of  ftj  values  was  obtained  by  using  a Fast 
Fourier  Transform  on  i reT^u^-1^l 

IT  1 J 

The  third  program  approximated  S (ft) 


, -T 

2e  r°° 


S2  (ft)  = f du  cos  (ftu)  f d£  e 

ttV2?  ° 0 


00  nr  «-^/2  r T<K?u) 


Here  tne  £ integration  was  performed  from  0 to  TMAX,  an  input  parameter.  The 
u integration  was  again  performed  using  a Fast  Fourier  Transform  to  get  a set 
of  answers  at  points  ftj  as  in  the  first  problem. 


where  JQ  is  the  Bessel  Function  of  zero  order.  Here  the  u ingegration  was 

T£2 

performed  from  zero  to  TMAX,  where  TMAX  = minimum  (UMAX,  , if  X^O) . Here 

X2 

UMAX  and  T£  are  input  parameters.  T £ is  generally  set  to  be  the  argument  of 
the  £th  zero  of  the  Bessel  Function.  The  value  of  £ is  also  read  in  and  the 
numerical  integration  starts  with  4£  + 1 equally  spaced  values  so  that  no 
oscillation  of  the  Bessel  Function  will  be  missed.  The  X integration  is  again 
performed  using  a Fast  Fourier  Transform  to  get  a set  of  answers  at  points 
ftj  as  in  the  first  problem. 

These  programs  lead  to  very  fast,  fairly  accurate  answers  (at  least  two 
significant  figures).  The  Initiator  indicated  that  these  results  were  helpful 
in  his  analysis  of  electro-magnetic  signals  in  the  ionosphere. 


VlglcodeA.  PteAwMi&Lon  ofi  fooUuoaJL  Angle. 

I onoiphenlc  Wave^omi 

Initiator:  Dr.  W.  Pfister 

Problem  No:  3028-6  Project  No:  7663 


A technique  for  measuring  drift  motion  and  turbulence  characteristics  of 
various  layers  of  the  ionosphere  consists  of  receiving  reflected  signals  from 
the  ionosphere  at  various  pulse  frequencies  in  the  MHz  ranges  of  closely 
spaced  antennas  and  digitally  recording  the  complex  amplitudes  of  the  echo 
pulses  in  the  time  domain.  This  is  accomplished  by  use  of  a Complex  Amplitude 
Multifrequency  Scanner  in  conjunction  with  a digital  ionosonde. 

The  time  domain  information  can  subsequently  be  spectral  analyzed  for 
power  content,  and  the  power  per  unit  area  can  be  presented  graphically  to 
represent  characteristic  ionospheric  fluctuations. 

The  purpose  of  this  problem  is  to  reconstruct  spectral  information  and 
to  present  a map-like  display  of  reflected  waveform  axial  projections  by  use 
of  a digicoder  printing  device. 
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A digicoder  display  consists  of  a continuous  printout  containing  4 lines 
of  data.  Each  line  consists  of  256  characters,  allowing  a 2-dimensional  map 
of  (256  x N)  characters  to  be  printed.  A unique  numerical  character  set, 
specific  to  the  digicoder,  allows  16  values  of  perceptive  weight  to  be  ob- 
served, thereby,  adding  the  dimension  of  color  density  (white  to  black)  to 
the  numeric  content  of  each  printed  character.  Thus,  the  points  in  a peak 
of  2 -dimensional  data  are  presented  as  a darker  area  relative  to  background 
information,  when  observing  the  total  map  of,  say,  (128  x 256)  points,  yet, 
the  specific  value  of  each  individual  point  is  retained  on  any  given  scale 
of  16-point  resolution. 

Input  to  the  computer  program  consists  of  Doppler  Frequency-related 
Power,  Coherence,  and  Phase  spectra  previously  calculated  from  the  Fourier 
Transform  of  the  time  domain  raw  data. 

Basically,  the  phase  relations  are  translated  into  angles  in  space,  as 
functions  of  the  recording  hardware  operating  frequencies  and  the  geometry 
of  the  receiving  antenna  orientations. 

The  locations  of  the  reflecting  regions  are  calculated  from  the  cross- 
spectral  phases  between  pairs  of  antennas,  at  those  Doppler  shifts  characterized 
by  peaks  in  the  Power  spectrum.  The  cross -spectral  phase  provides  an  estimate 
of  the  direction  of  the  reflecting  region.  A band  of  width  A<p  of  the  reflect- 
ing region  is  related  to  the  given  spectral  coherency  (COH)  according  to: 

A4>  = ± tt(1->/55h) 


The  intersection  of  2 bands  from  2 antenna  pairs  provides  a parallelogram 

area : 


A = — — (1-v/C^D  (1-v/EofT) 

Via2  1 2 

where  X = the  transmitter  operating  wavelength 
and  a = the  length  of  the  antenna  pair  baseline 
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The  power  density  per  unit  area  is  defined  as  W/A,  where  W in  this  program 
is  the  average  power  from  3 receiving  antennas  at  a particular  Doppler 
frequency. 

From  the  geometry  of  the  3 receiving  antennas,  three  such  parallelograms 
can  be  constructed,  and  from  each  parallelogram,  trapezoidal  areas  can  be 
projected  onto  X and  Y axes  oriented  in  the  N-S  and  E-W  directions.  The  sum 
of  the  heights  of  3 projected  trapezoidal  areas  on  a given  axis  is  proportional 
to  the  spectral  density,  W,  and  is  indicative  of  the  strength  of  the  power 
content  of  the  given  region.  The  angular  location  of  the  reflecting  region 
(see  mathematics  below)  is  given  from  the  range  over  which  each  summed 
trapezoidal  area  extends  on  each  projection  axis.  From  the  sums  of  the 
trapezoidal  projections,  then,  two  intensity  modulated  lines  can  be  displayed 
on  the  digicoder,  each  as  a function  of  Doppler  frequency  and  arrival  angle. 

It  should  be  noted  that  the  program  contains  test  criteria  to  avoid 
phase  ambiguities  in  the  reflected  signals.  It  was  desired  that  the  range 
of  the  display  correspond  to  the  principal  range  of  phase  between  2 antennas, 

-it  to  it  radians.  The  test  criteria  are  dependent  upon  the  antenna  pair  chosen 
for  observation  and  the  projection  axis. 

Each  digicoder  output  map  is  specific  to  a given  time  of  recording,  which 
is  contained  in  a preface  (identification)  area  preceding  each  powpr  density 
display.  In  the  digicoder  printout,  128  points  of  Doppler  frequency  are 
presented  along  the  vertical  (length)  of  the  display  and  256  points  (128  per 
N-S  and  E-W  projection  axis)  along  the  horizontal  (width).  In  addition, 
markers  are  presented  every  16th  character  along  the  horizontal  in  the  preface 
area  to  assist  in  notation  of  scale.  In  the  same  manner,  preface  identifica- 
tion characters  are  presented  every  16^  line  along  the  left  side  of  the 
vertical  to  denote  the  scale  of  the  Doppler  frequency. 

Overall,  for  one  subcase  of  time,  four  such  (128  x 256)  point  maps  are 
generated,  each  corresponding  to  a given  height  range  gate  on  the  recorder. 


B 
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Hie  program  reads  input  from  a packed  magnetic  tape  generated  at  a prior 
stage  of  production  processing.  This  program  then  generates  digicoder 
information  on  packed  magnetic  tape,  to  be  used  at  a later  stage  of  data 
processing  as  digicoder  input. 

Mathematics 

1 Location  of  reflecting  areas  on  X-  and  Y-axis: 


Antenna  Pair 

X Axis  Location 

Y Axis  Location 

12  , 23 

- <P  - <t>  * 2ft/ 100 

— (-  <t> 

+ <|>  ) 

12  23 

n/3 

2 3 

23  , 31 

<t>  + 2ft/ 100 

— (2  <P 

+ * ) 

3 1 

y/J  23 

31 

31  , 12 

♦ + 28/100 

— (-  <t> 

+ 2<J> 

3 1 

y/J  31 

12 

where  $12,  <p23  and  <p31  represent  cross-spectral  phase  from  antenna  pairs 
12,  23,  31. 

and  ft  = Doppler  frequency  under  investigation. 

2 Maximum  Spectral  Intensity  projected  onto  trapezoidal  areas  along  the 


E-W  (Y)  axis. 


Area  Reference  No. 


Intensit 


1 


V3 


vT 


A<|>  + 2A<j>  + 1^1  A<|>  - 2A<p  I 


31 


12 


31 


12 


where  Q - | ^ 

W = average  power  at  a particular  Doppler  frequency 
X = operating  wavelength  (meters) 
a = antenna  baseline  (meters) 

3 Maximum  Spectral  Intensity  projected  onto  trapezoidal  areas  along  the 
N-S  (X)  axis. 

Area  Reference  No.  Intensity 
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VoppteA  Ski^t/lonoAphe/uc  Rejection 


Initiator:  Dr.  W.  Pfister 

Problem  No:  3028-7  Project  No:  8658 


This  program  is  part  of  a data  processing  project  called  DAASM  (Doppler 
Angle  of  Arrival  Spectral  Measurement),  for  detailed  study  of  the  ionosphere. 
This  particular  effort  involves  the  estimation  of  Doppler  shift  and  the 
combining  of  three  Fourier  transforms.  Signals  generated  in  this  case  by  a 
moving  aircraft  are  reflected  by  the  ionosphere  and  measured  on  the  ground. 

The  results  are  recorded  on  magnetic  tape  as  a time  series  of  digitized 
complex  numbers.  The  computer  program  reads  the  magnetic  tape,  isolates 
samples  about  one  minute  long,  calculates  three  overlapping  Fourier  transforms 
from  each  sample,  and  attempts  to  determine  the  Doppler  shift  of  the  strongest 
frequencies  by  looking  for  a consistent  phase  shift  among  the  three  transforms. 
Then  the  three  transforms  are  combined,  taking  into  account  the  Doppler  shift. 
The  results  are  written  on  two  magnetic  tapes,  one  in  a format  suitable  for 
the  Digicoder  plotter,  the  other  for  input  to  other  programs. 

The  raw  data  input  tape  is  pre-positioned  at  the  desired  starting  case. 

An  input  card  is  read  specifying  the  number  of  consecutive  cases  to  be  pro- 
cessed. Each  case  occupies  seven  physical  records  on  the  input  tape. 


[ 

i 

t 

: 
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Each  case  is  then  handled  as  follows: 


The  first  record  is  skipped.  Prefaces  are  extracted  from  the  remaining 
records  and  printed.  Record  numbers  are  checked.  If  normal  sequence  is 
broken,  the  program  starts  again  from  the  current  record.  Data  are  packed 
into  the  array  BUF.  Then,  taking  one  frequency  and  one  antenna  at  a time, 
we  obtain  three  Fourier  transforms  from  the  three  interlaced  complex  data 
sets.  (The  three  interlaced  sets  were  all  sampled  the  same  way;  sets  2 and  3 
are  just  B = 6/25  of  a step  later  than  sets  1 and  2 respectively,  where  one 
step  is  the  time  between  consecutive  samples  in  one  set,  .25  seconds). 
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Each  Fourier  transform  contains  256  complex  points  gm,  m = 1,2,..., 256, 
each  point  representing  the  summed  contributions  of  a whole  family  of  aliased 
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frequencies  (m-1  ♦ 256  n)/64  Hertz,  n = 0,  ±1,  ±2,...  The  problem  is  to 
pick  the  n corresponding  to  the  "best"  (usually  strongest)  frequency  in  the 
original  sampled  signal.  For  this  purpose  we  assume  that  a singly  value 
of  n accounts  for  all  the  aliasing  in  the  neighborhood  of  the  best  frequency. 
We  divide  the  aliased  frequency  range  into  four  equal  segments  and  analyze 
each  segment  to  see  which  best  fits  the  model. 

A displacement  of  At  seconds  in  the  time  series  sample  causes  phase 
shifts  proportional  to  frequency  in  the  Fourier  transform.  If  we  assume  as 
above,  that  entry  gm  in  the  transform  is  all  due  to  one  frequency,  namely 
(m-1  + 256  n)/64  Hertz  (n  is  a fixed  but  unknown  integer),  then  a At  displace- 
ment will  cause  a phase  difference  of  At  • (m-1  + 256  n)/64  cycles.  And,  if 
the  same  n applies  to  a sequence  of  gm's,  the  phase  differences  between  two 
interlaced  transforms  will  be  a linear  function  of  m. 

Let  jgnk,  m=l,2, . . .256 } be  the  transform  of  the  k**1  interlaced  time 
series,  k=l,2,3.  The  following  expression  uniquely  defines  p,,^: 


gmk  = HgnJI  exP  t271^’0  - Pmk  < 1 


We  have  At  = .250  seconds,  so,  if  our  model  were  correct  we  would  have: 


oyk.i  - p„k>  „ - ce[o.-i)/2S6«j) 

modi  modi 


or 


def 


"mk 


<P».k*l  - Pmk  - 8<»-W256>»odl  ' <6">modl 


In  other  words,  the  c,^  should  all  be  about  the  same  and  will  tell  us 
what  n is. 

In  each  segment  of  the  transform,  (m  = 1-64  or  65-128  or  129-192  or 
193-256) , we  calculate  the  weighted  average  A of  unit  vectors  with  argument 
I’TCmk  on  the  complex  plane.  The  weighting  is  with  respect  to  power. 


11 


A = ( E E 
k=l,2  m 


lgmk  «m,k+lH  exP  (2lricmk))/(  2 2 M*mk  ^.k+l^ 

* k*l,2  m * 


= ( 2 2 IIC  gm  k+l  exp  [-2iri6(n-l)/256]||)/(  E E ||g  g k ||) 

k-1,2  m nnK  m’K  1 k-1,2  m 1 


We  also  calculate  in  each  segment  an  estimate  of  the  scatter  of  the  c^. 


a = (1-AA  ) 


‘/2 


We  pick  the  segment  with  the  smallest  scatter.  The  corresponding  value  of  A 
yields  our  best  estimate  of  n;  we  pick  an  n minimizing: 


1 1 A/ 1 1 A 1 1 - exp  (2Tii0n)  | j . 

Since  3 = 6/25  and  n is  an  integer,  it  suffices  to  consider  n = 1,2,..., 25. 

With  n and  the  best  segment  well  determined,  the  three  interlaced 
transforms  may  be  combined  in  a meaningful  way.  Let  m^  be  the  final  subscript 
of  the  best  segment,  m^  = 64,  128,  192  or  256.  The  new,  combined  transform 
will  cover  the  frequency  range: 

(n^  - 128  ♦ 256  n)/64  Hertz 


to 


(mjj  + 127  ♦ 256  n)/64  Hertz. 


We  arbitrarily  assume  that  the  phase  shift  is  linear  over  that  entire  range 
and  that  the  g^'s  are  all  aliased  from  that  range.  The  phase  of  each  gmi  and 
gm3  is  modified  to  align  with  g,,^,  and  the  three  are  averaged.  The  results 
are  plotted  and  the  digicoder  with  the  (mb  - 1 ♦ 256  n)/64  Hertz  line  is  in 
the  middle  (i.e.,  the  128^  point). 


* 
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This  method  breaks  down  for  cases  in  which  the  original  signal  contains 
two  good  frequencies  which  cannot  without  aliasing  be  included  in  a single 
2S6  point  range.  Only  one  of  the  peaks  will  be  smoothed  correctly. 

NOTE:  The  Fourier  transform  subroutine  used  by  this  program  delivers 

the  usual  backward  transform.  The  write-up  above,  however,  is 
from  the  point  of  view  of  the  usual  forward  transform,  which 
gives  a larger  phase  to  a signal  sampled  later.  Thus,  the 
program's  G arrays  are  the  conjugates  of  the  g's  in  the  write- 
up, the  original  signal  also  having  been  conjugated. 

In  other  words,  given  a time  series  f , ...,f  , subroutine  FORER 

1 25o 

calculates  gj,...,g2S6 


®m+l 


exp  (imj  2tt/256) 


of 


exp  (-imj  2tt/256)  . 


The  computer  program  outputs  a heading  including  the  input  value  of  the 
maximum  number  of  cases  to  process. 

Then,  for  records  2-7  of  each  case,  it  prints  the  prefaces.  For  each  of 
the  two  frequencies,  first  antenna  only,  these  are  followed  by  the  segment 
number  of  the  best  segment,  the  best  value  of  n in  that  segment,  the  point 
number  in  the  original  Fourier  transform  which  will  appear  on  the  extreme  right 
when  re-arranged  for  the  Digicoder  plot,  and  the  number  of  Hz,  module  100, 
at  the  center  of  the  Digicoder  plot.  Then  come  the  gn  averages  and  associated 
a estimates  for  each  of  the  four  segments. 
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Input  Tape  Format 

Each  case  occupies  seven  physical  records.  Each  physical  record  is 
343-344  words,  comprising  3430  to  3440  6-bit  bytes.  Before  the  input  can  be 
interpreted,  it  must  be  aligned  properly.  Let  k be  the  minimum  integer  such 
that  250  < k < 323  and  such  that  the  k^  byte  of  the  record  is  the  beginning 
of  a 288  x 11  array  of  bytes  with  all  zeros  in  the  first  18  rows.  If  there 
is  no  such  k,  set  k to  271.  Slide  the  entire  record  so  that  the  kth  byte  will 
be  in  the  271st  position.  If  k < 271  fill  the  initial  part  of  the  buffer  with 
zeros.  Then  the  following  format  applies. 

Consider  the  bytes  to  be  arranged  in  standard  Fortran  order  in  an  array 
dimensioned  (3,  2,  3,  4,  48).  (The  last  few  positions  of  this  array  will  be 
unoccupied).  Then,  for  £ = 1,  2,  3,  the  (i,  j,  k,  £,  m)-th  byte  of  the 
array  is  the  ith  byte  of  the  mth  sample  from  the  £ih  interlaced  set  from  the 
kth  antenna  at  the  jth  carrier  frequency.  The  three  bytes  (i  = 1,2,3)  in 
each  sample  are: 

byte  1:  |Re(z)| 

byte  2:  1 in  the  3r^  bit  (10  octal)  if  and  only  if 

Re(z)  >_  0,  plus  1 in  the  4^  bit  (4  octal) 
if  and  only  if  Im(z)  0. 

byte  3:  | Im(z) | 

For  £ = 4,  m > 3,  there  is  no  data. 

For  £ = 4,  m = 1,2,3,  there  are  18  bytes  available  for  each  of  three 
"prefaces". 


1 - bytes 

1-3 

Code  digits  6,  7,  3 

bytes 

4-6 

Three  decimal  digits  of  day  of  year 

J 

bytes 

7-12 

Six  decimal  digits  of  hours,  minutes, 
and  seconds. 

* 

bytes 

13-14 

Code 

•j 
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hi  ■ 


ri 


' i 

I 

I 

r 

i 


I 

! 


i 

! 


I 


i 

1 


Two  digits  containing  double  the 
record  number  within  this  case. 

Code 

Three  decimal  digits  giving  frequency 
setting  of  Is*  antenna,  1st  frequency. 

Same  as  1-3  but  for  2nd  frequency. 

Same  as  1-6  but  for  2nd  antenna. 

Same  as  1-6  but  for  3rd  antenna. 

Same  as  Preface  2,  except  3 digits 

give  range  setting  instead  of  frequency 
setting. 

There  is  an  end-of-file  after  the  last  record  of  the  last  case. 

If  the  first  byte  of  a frequency  is  0,  1,  6,  7,  10,  11,  12,  or  13,  the 
corresponding  raw  sample  must  be  conjugated  to  compensate  for  a hardware 
difference  in  the  way  the  data  were  collected. 

Output  Digicoder  Tape  Format 

The  output  tape  is  formatted  for  the  Digicoder,  a machine  which  produces 
shaded  plots. 

The  program  writes  two  files,  one  for  each  frequency.  Each  file  has  a 
header  record,  repeated  3 times,  then  one  data  record  for  each  case  on  the 
input  tape.  The  header  record  contains  preface  information  in  various  standard 
arrangements. 

Each  record  has  245  words,  or  2450  6-bit  bytes.  The  last  two  bytes  are 
zero.  The  other  2448  are  a 408  x 6 array,  plotted  in  six  lines.  The  first 
24  bytes  of  each  line  are  a preface,  the  same  for  all  six  lines,  as  follows. 


Preface  1 - bytes  15-16 
(Cont . ) 

bytes  17-18 

Preface  2 - bytes  1-3 

bytes  4-6 
bytes  7-12 
bytes  13-18 

Preface  3 - bytes  1-18 
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Byte  2: 

B j 

Byte  3: 

f 


50  (octal)  if  and  only  if  vj  0 plus  4 
(octal)  if  and  only  if  V2  >_  0. 


Line  1 contains  the  real  parts  of  the  combined  antenna  1 transform; 
line  2,  the  imaginary  parts.  Lines  3-6  contain  the  transforms  from  antennas 
2 and  3,  with  the  same  Doppler  shift. 


Combined  Aircraft  Data  Transform  Tape  Format 


This  tape  contains  one  file  consisting  of  two  logical  records  for  each 
case  processed.  There  is  one  record  for  each  carrier  frequency.  Each  record 
contains  1591  words  as  follows. 


WORDS 


1-18 

LPX  array  of  18  digits  from  Preface  1 
of  the  first  record  used  in  this  case. 
Integers. 

19-36 

FREQ  array  of  18  digits  from  Preface  2 
as  above . Integers . 

37-54 

RNG  array  of  18  digits  from  Preface  3 
as  above . Integers . 

55 

Hz  the  number  of  Hertz  at  the  129^ 
member  of  the  following  transform.  Given 
module  100;  -6  appears  as  94. 

56-1591 

Complex  array  H in  normal  Fortran  order, 
dimensioned  256  x 3.  H(I,J)  = the  combini 
transform  entry  for  the  Jth  antenna  at  Hz 
+ (I-129)/64  Hertz. 

Sky  PolaJiLzaution 
Initiator:  F.  Volz 

Problem  No:  3034-5  Project  No:  7621 


A series  of  four  programs  were  written  to  aid  in  the  reduction  of 
observations  >f  neutral  points  of  sky  polarization  in  the  clear  sky  during 
twilight.  The  Initiator  has  data  gathered  over  a three  year  period,  showing 
evidence  of  changes  in  relation  to  variations  of  stratospheric  aerosol 
content . 


The  main  portion  of  the  effort  was  directed  towards  presentation  of 
the  data  in  a variety  of  printed,  punched,  and  plotted  formats.  Certain 
corrections  and  interpolations  were  performed.  In  addition,  it  was  necessary 
to  write  a subroutine  calculating  solar  elevation  as  a function  of  date, 
Universal  Time,  and  station  latitude  and  longitude. 
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Vamping  Functions  In  Infrared  Absorption 

Initiator:  Dr.  Bernard  Bendow 

Problem  No:  3055-3  Project  No:  3326 


This  problem  is  one  of  a set  of  problems  directed  toward  the  prediction 
of  the  infrared  absorption  intensity  and  line  shape  as  a function  of  frequency 
for  various  models  of  ionic  crystals  and  semiconductors.  In  particular,  this 
problem  continues  the  work  performed  under  Problem  No.  3055-1  for  Air  Force 
contract  F19628-70-C-0120. 

This  problem  required  the  calculation  of 


!?,  R1  ,R  n_1  30  D=0 

3 

, yi  refers  to  a sum  over  all  the  latice  sites  of  a unitary  cubic 

R.R1,^ 


1 rBF 


where 


crystal,  and  where  distances  are  measured  from  one  of  the  latice  sites. 


Here : 


pn+l(u)  = £°  d“1  pn(“"a,1) 


p (to)  = C (-!■*>)  + C (R1  -R-]t  ,o))  - C (R^-R1,^  - C (It-R  ,w) 


n 3 


12  3 


F = exp  (-  ^5® ) (E-3) 
E2 


where , 
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where 


A = a ( | R| ) 2 + a (|R|)2;  B = R-Rl;  E = (4a  a -D2)1 
1 2 1 12 


R = kn  ; R1  = R*+n  ; a = i(C  (0,0)  ♦ C (0,0)  - C (0,0)) 
' 1 2 11  22  12 


+ Ro2;  a = -(C  (0,0)  + C (0,0)  - C (0,0))  + Ro2 

2 2 11  22  12 


where,  the  vector  n is  an  input  parameter. 

An  explanation  of  the  general  theory  leading  to  equations  of  this  form 
was  also  presented  in  the  final  report  for  the  above  mentioned  contract. 

For  this  problem  two  forms  of  C^j  were  tried.  These  two  forms  are 


C.  . (£,w)  » — n,  4 m ” 1 
ij  ’ ^2  b 2 s 


(w-w  y)  . 

sin  [Ko  THaJy  ,R|] 

jtojij 


eiek  C“  1 

,/nr TT^  * 11  £orY“1i“iY 

V J X w(l-u)  ) e -1 


Cij(lf,a))  = e^  C (It.  -w)  for  -y  < u £ 


Ci^(R  ,to)  = 0,  otherwise 


where 


ei  \m1*m2  * 


(f"m  m +m  m +m 

— i — r = _i — L R = -A — - 
m ♦m  i in  2 m 

12  2 i 


k- 


and  y,  nb.  ms»  Ko»  > “j * m2»  0 are  inPut  parameters 


i&K-OBOiHeBnaBi 


and 


C.  . (R,u>)  = -n,3/z  m _1  ~ .inlKo-Kll -RJ I i (-J — +1)  for  0 < <*>  < y 

" Y d s Vs A 2C1-'° " “ 


l.l 


Ci . (R,w)  = (R,  -w)  for  -y  < u < 0 


(R,gj)  = 0,  otherwise 


where 


K = 1 - ■yi-  y for  0 U)  < 1 K = 1-DELTA  if  w=l; 


Ko  = (6ir2n  )x^3  and  y,  n,  , m , $,  n , DELTA,  m , m 

C DSC  2 2 


are  input  parameters. 

To  perform  the  calculations,  a program  was  written  which  did  the  algebraic 
manipulation  needed  to  derive  the  derivatives  of  F.  This  program  generated 
card  output  for  insertion  in  other  programs,  which  do  the  actual  calculations. 
Derivatives  up  to  the  twelfth  order  were  obtained.  This  order  of  differentia- 
tion was  found  to  be  appropriate  in  evaluating  I(w). 

The  computer  program  used  the  trapezoidal  rule  to  approximate  the 
integration  needed  to  generate  pn  from  pi  using  the  recursion  relationship 
Pn+l(w)  = fZ  dw1pnCu-w1)  PjCu1)  and  the  knowledge  that  Pn(t)  = 0 for  |t|  > ny. 
The  program  calculated  and  plotted  results  of  I(ui)  versus  u for  all  the  values 
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of  the  sum,  such  that,  | R | 2 ♦ f R1 1 2 < IR  (an  input  parameter)  and  | < M, 
where  M equals  1,2,3,...,  ISRRT  (another  input  parameter). 

A naZy-i-U  o&  the.  TempeAatu/ie.  Gradient  in  tht  ktmo^phzne. 

Initiator:  Mr.  A.  Cole 

Problem  No:  3065-1  Project  No:  8624 

The  purpose  of  this  problem  is  to  provide  an  objective  method  for 
obtaining  the  horizontal  temperature  gradient  between  the  North  Pole  and  the 
equator  from  observed  data  at  levels  between  30  and  90  km.  In  essence,  this 
method  smooths  the  observed  temperatures  and  provides  interpolated  values 
for  15°  intervals  of  latitude.  The  interpolated  temperatures  were  to  be  used 
with  appropriate  atmospheric  pressures  to  construct  U.S.  Reference  Atmospheres 
representing  mean  monthly  conditions  at  15,  30,  45,  60,  75  and  90°  N.  latitude 

For  a given  month,  the  Initiator  provided  13  sets  of  data.  Each  set 
representing  measurements  of  temperature  versus  latitude  for  a fixed  altitude. 
The  method  of  least  squares  is  used  to  determine  the  best  fit  for  each  data 
set  by  polynomials  of  order  1 thru  6.  The  results  of  this  processing  were 
presented  in  a tabulated  form  that  facilitates  comparison  between  the  order 
of  polynomials  fitted  to  a particular  data  set.  This  comparison  is  essential 
to  the  Initiator  to  ascertain  his  preference  that  for  a specific  altitude 
and  month,  a certain  polynomial  best  describes  the  dependence  of  temperature 
on  latitude.  This  effort  utilized  the  least  squares  program  of  the  Analysis 
and  Simulation  Branch  Computer  Library. 
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Expomntiat  CuAv&  F-btUng  faoK  UV  Abio^ption  Data 


* 


Initiator:  Mr.  L.  Weeks 

I 

Problem  No:  3069  Project  No:  6690 

I 

i 

I 

The  purpose  of  this  problem  is  to  analyze  ultra-violet  absorption  data  *j 

obtained  in  the  110  to  220  km.  region  from  rocket-borne  photometers.  The  V 

results  of  this  analysis  are  in  turn  utilized  in  the  determination  of  the 
molecular  oxygen  number  density  of  this  region. 

i 

The  expected  behavior  of  ultra-violet  absorption  data  (I(h^)  is  that  it 
undergoes  an  exponential  variation  with  respect  to  altitude.  This  behavior 
requires  the  determination  of  the  'best'  values  for  the  parameters  IQ,  hQ,  Ht 
and  H2,  such  that,  the  measured  data  I(hj)  can  be  represented  by  the  expression 

I 

I(hi)  = IQ  exp  - j.5  exp  - 0^-h^/H^  + .5  exp  - O^-hp/H^ 

where,  hj  is  the  independent  variable  for  altitude  and  'best'  values  are  in  j 

accordance  with  the  criteria  of  least  squares,  i.e.,  minimize 

l 

i 

<J>  = E (I(h^)  - Expression)2  1 

i 

Using  standard  iteration  techniques,  convergence  was  not  obtained  for  the 
curve  fitting  in  the  prescribed  form.  Also,  convergence  was  not  obtained  when 
a logarithmic  transformation  was  applied  to  both  the  data  and  the  prescribed 
form  of  the  curve  fit. 

The  'best'  estimates  for  the  parameters  were  successfully  obtained  when 
utilizing  the  Levenberg  Method  and  simultaneously  constraining  the  parameters 
H and  H . A full  description  of  this  method  and  its  application  are  presented 

1 2 i 

in  the  Final  Report,  AFCRL-TR-73-0433  mentioned  in  the  introduction  of  this  .! 

report . 

The  results  of  this  processing  were  presented  at  the  AGU  conference  in  * 

San  Francisco  in  December,  1973. 


Tracking  and  AnaZyAiA  Variable.  Low  Frequency  Uavt^onm 

Initiator:  Capt.  McLain 

Problem  No:  3070-1  Project  No:  4603 


Program  PROGS  was  written  in  support  of  a study  of  Variable  Low  Frequency 
Waveform  Tracking  and  Analysis.  The  purpose  of  this  segment  of  the  study  is 
to  calculate  and  plot  Ionospheric  Reflection  Coefficients  and  Group  Heights 
from  VLF  ionospheric  waveform  reflections  recorded  as  a function  of  time. 

Logically,  the  computer  program  proceeds  in  the  following  manner: 

Sequences  of  three  pulses,  identified  as  groundwave,  skywave  and  rotated 
skywave,  are  unpacked  from  magnetic  tape  and  Fast  Fourier  analyzed  for  frequency 
content.  Each  Fast  Fourier  spectrum  is  fitted,  according  to  a third  order 
polynomial,  to  a given  frequency  interval,  and  the  resultant  amplitude  and 
phase  values  of  two  predetermined  frequencies  are  selected  for  use  in  computa- 
tion of  reflection  heights  anti  coefficients.  An  amount  of  averaging  of  fitted 
Fourier  data  takes  place  in  the  time  domain  according  to  an  averaging  time 
span  specified  by  the  user. 

J 

For  each  given  time  period  to  be  analyzed,  a series  of  10  pairs  of  pen- 
and-ink  plots  are  generated,  consisting  of  8 pairs  of  families  of  phase  heights 
and  coefficients,  and  2 pairs  of  group  heights  and  coefficients,  all  presented 
as  functions  of  time  (in  hours).  All  phase  heights  are  plotted  over  6 inches 
on  a linear  scale  of  0-120  km.,  and  reflection  coefficients  are  plotted  over 
4 inches  on  a logarithmic  scale  encompassing  the  two  decades  between  0.01  and 
1.0.  Each  plot  is  identified  according  to  day,  year,  time,  frequency  of 
analysis,  and  rotated  or  normal  skywave  reference  data. 

The  user  has  control  over  several  parameters  within  the  program  via  input 
control  cards.  Included  among  these  are  the  number  of  time  series  to  be  analyzed,  I 
the  spacing  at  which  the  Fast  Fourier  spectra  are  to  be  fitted,  a 30  character 
plot  label  array,  the  amount  of  time  over  which  fitted  Fourier  data  averaging 
is  to  take  place,  the  starting  and  ending  times  over  which  analysis  is  to  take 
place,  and  the  two  frequencies  selected  from  the  fitted  Fourier  spectra  at 
which  analysis  is  to  take  place. 


Mathematics 


A quantity  of  interest  is  the  phase  difference  between  the  ground  and 
skywaves  (PHSG) . 

The  mirror  height  is  calculated  from  PHSG  by  the  following  formula: 


HTMIR  = j-\[DG+(PHSGC*N*T)  * .29979]^  - 


DG  * distance  between  transmitter  and  receiver  (KM) 

PHSGC  = PHSG-PHM+PHGC  (usee) 

PHSG  ■ (skywave  phase  data  - groundwave  phase  data)  - (Skywave  phase 
calibration  - groundwave  phase  calibration). 

0<PHSG<T 

PHM  = Phase  shift  of  the  skywave  at  the  ionosphere  (ysec) 

PHM  _ or  T/2  for  the  normal  component 

(T/4  or  3T/4  for  the  rotated  component 

PHGC  = Correction  to  groundwave  phase  at  frequency  F;  due  to  the  finite 
conductivity  of  the  ground  (ysec) 

N = Number  of  cycles 

0<N<*  (Upper  boundary,  N*,  is  established  when  HTMIR  exceeds  120  KM), 

T = 1000. /F  (T  in  ysecs,  F in  kHz) 

.29970  = Velocity  of  Light  (KM/ysec) 

Reflection  Coefficients  - The  reflection  coefficient  is  a ratio  of  the  skywave 
to  groundwave  signal  strength  with  certain  corrections.  The  reflection 
coefficient  is  calculated  from: 


RCOEF  = 


2RAC*AMGC* 


,2tt*F*AL*cos0.1  2tt*F*AL  % 

cos  ( 5—)  - cos  (- «■) 

9.83571x10  $.83571x10* 


f in  r nil  . 
COS  ( -.) 

9.83571x10^ 


where 


ARSG  = component  signal  amplitude 

Groundwave  signal  amplitude 


or 


in  (Skywave  amplitude  (dB)  - Groundwave  amp  (dB) 
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if  the  signal  amplitudes  are  given  in  decibels 


RAC  = Receiver  antenna  response  amplitude  correction 

<1.  For  the  normal  component 
( CosQ  for  the  rotated  component 

AMGC  ■ Correction  to  the  groundwave  amplitude  at  frequency  F due 
to  the  finite  conductivity  of  the  ground 

AL  * Antenna  length  (KFT) 

0 ■ Skywave  angle  of  incidence 

0 = Tanl^2HTHlT^ 

9.83571xl02  = Velocity  of  Light  in  KFT/psec 
F * Frequency  (kHz) 

Reflection  heights  and  coefficients  are  calculated  for  both  mirror  phase 
shifts  of  both  skywave  components  at  each  frequency.  These  values  are  computed 
for  each  data  record  interval  and  plotted  vs.  time  of  day. 
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Symbolic.  Solution  o&  a Linear  Equation 

Initiator:  David  Anthony 

Problem  No:  3074-1  Project  No:  7600 

This  problem  involved  solving  for  the  C in  the  following  equations: 

NVA 

aIJ  C(2J.O)  Kj  for  I “ 1 NVA 

where  ajj  and  Kj  are  polynomials  in  two  variables  F and  Ml  and  the  solution 
is  to  be  in  terms  of  F and  Ml. 

The  method  used  to  solve  these  equations  is  to  reset  the  aTT  and  KT  to 

1J  1 

new  polynomials  as  follows: 

Reset  for  1=1  then  I = 2 I * NVA 

ajK  is  reset  to  an  aj|C  - a.j  aIK 

j t I 

K * I 

K.  is  reset  to  aTT  K.  - a.T  K_ 

J II  j jl  I 

j t I 

a^ j is  reset  to  0. 

j * I 

The  resulting  values  of  C(2I,0)  are  then 


C(2I,0)  KI  = 

Here  Kj  and  a^  are  polynomials  in  F and  Ml. 

These  results  were  then  expanded  in  powers  of  F so 


where  each  term  is  a polynomial  in  Ml. 


'(21,0)  ‘ a 


II 


F=0 


NK-1 

l (( 

n=l 


3h(~) 

aII 

3h  F 
TT 


f=0 


26 


V^ociatLon  PxeA&uAjt  of  Silicon  CoKbi. te. 

Initiator:  Ur.  J.  Smiltens 

Problem  No:  5077-01  Project  No:  5620 

I 

The  equilibrium  pressure  for  silicon,  carbite,  grafite  and  vapor  has  been 
determined  by  two  groups  of  investigators.  In  1930  in  Germany,  Grieger 
determined  it  in  the  range  from  760  mm.  mercury  to  1 mm.  mercury.  Grievenson 
and  Alcook  determined  it  at  lower  temperature  in  the  micropressure  range. 
According  to  classical  Physical  Chemistry,  the  logarithm  of  vapor  pressure 
when  plotted  vs.  reciprocal  absolute  temperature  should  give  a straight  line. 
Both  investigators  have  obtained  the  straight  line.  However,  Grieger' s straight 
line  has  a larger  negative  slope  than  Grievenson' s and  Alcook' s.  A more 
careful  consideration  and  recent  experimental  evidence  indicate  that  Grieger 's 
data  should  depart  upward  from  the  Grievenson* s and  Alcook' s straight  line. 

The  special  condition  the  parabola  should  be  tangent  to  the  straight  line. 

A functional  description  of  the  required  mathematical  analysis  and  computer 
programming  is  the  following: 

Calculate  the  least  squares  polynomial  through  the  submitted  Grievenson 
and  Alcook' s data. 

Calculate  the  least  square  parabola  through  the  submitted  Grieger' s 
data. 

Solve  a system  of  the  submitted  two  equations  with  C and  XQ  unknowns. 

Find  the  real  root  of  the  resulting  5th  degree  polynomial  which 
satisfies  the  special  condition:  the  parabola  should  be  tangent  to 
the  straight  line. 

» 

| 

The  solution  of  the  system  of  the  two  equations  submitted  with  C and  xQ 
unknowns  and  with  n representing  the  summing  of  the  data: 

2 £ [y  - A - Bx  - C (x* -2x  x ♦x*) ] (x*-2x  x +x*)  = 0 
4 n n o onnJ  o onn' 

2 l [y  - A - Bx  - C(x*-2x  x +x*)]  (2x  -2xJ  = 0 
w n n o o n n'J  ' o n 


y 
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omtt 


■u 


is  a S**1  degree  polynomial  a,  + a x ♦ a x2  a x3  + a x4  ♦ a x5  * 0 with 

o 1 2 3 4 s 

a = A6B5  ♦ A7B4 

0 

a » A4B4  - ASB5  ♦ A6B2  - A7B3 

1 

a * A7B1  - ASB2  - A4B3  + A3B5  + A2B4 

2 

a = A7B0  + A4B1  + A3B2  - A2B3  - A1B5 

3 

a = A2B1  * A4B0  * A0B5  - A0B2 
% 

a = -A0B2 
s 

where 

AO  = -n 
A1  = 4 Ex 

n 

A2  = Ey  - nA  - BEx 
7n  n 

A3  = 6Ex2 

A4  = 2(AExn  - Exnyv  ♦ BEx2) 

A5  = 4Ex2 
A6  = Ex4 

A7  = Ex2y  - AEx2  - BEx2 
n n n n 

BO  = AO 
B1  = 3Exn 

B2  = A2 

B3  = 3Ex2 
n 

B4  - Ex2 
B5  = A4 

NOTE:  AIBJ  = (AI)  (BJ) 
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R o&lional  Constant*  diatomic  Molecule* 

Initiator:  Dr.  Daniel  Katayama 

Problem  No:  3080-1  Project  No:  8627 


From  the  analysis  of  the  band  spectra  of  diatomic  molecules  in  the  vacuum 
ultra-violet,  the  rotational  constants  of  these  molecules  can  be  determined 
when  they  are  in  highly  excited  states.  These  constants  are  determined  by 
solving  a system  of  simultaneous  equations,  using  matrix  and  least  squares 
techniques.  Also  determined  are  the  variance  6Z  , and  the  variance-covariance 
matrix  (V)  described  below  in  the  following  equations: 


xi  - w1)  - A"2 


X"  = JV(JV.l)  - A"2 
where  J',  J",  A*  and  A"  are  data  input  numbers. 
The  system  of  simultaneous  equations: 


y,  = a + a (x‘)  - a (x') 

1 O li  2 1 


a (xV)  ♦ a (x"V 

3 1 4 1 


y = a + a (x*) 
2 o 12 


a (x')  - a (x")  <►  a fx"): 

2 2 3 2 4 2 


y„  = a + a (x'l 
7n  o i n' 


* au(xn)2 


Or,  in  matrix  form 


Corutnuctcon  otf  "RexuonabZe." , Analytic 
CuAveA  0(5  Wind  I/a.  Height  6*. om  Oit>cAe£e  Data. 

Initiator:  Rene  V.  Cormier 

Problem  No:  4026-6  Project  No:  8624 


Research  in  support  of  Air  Force  paradrop  operations,  investigating 
the  characteristics  of  vertically  integrated  boundary  layer  winds  established 
a need  for  the  construction  of  "reasonable"  analytic  functions  (curves)  of 
wind  and  temperature  versus  height  from  discrete  data.  These  data  were  obtained 
at  non-regularly  spaced  levels  ranging  in  height  from  the  surface  to  1500  ft. 
and  in  number  from  7 to  12  depending  on  the  tower.  The  curves  had  to  be  capable 
of  integration  and  differentiation,  and  because  of  the  number  involved  (over 
30,000)  had  to  be  generated  by  computer  without  human  intervention. 

A series  of  computer  programs  were  written  in  two  (2)  groups  and  were 
documented  under  the  name  C0RM  for  the  SUYA  Computer  Center  Library. 

GROUP  I GROUP  II 


1. 

CORMN 

5. 

CORMP 

2. 

INTGRLN 

6. 

INTGRLP 

3. 

HERMITN 

7. 

HERMITP 

4. 

C0RM4 

Programs  CORMN,  CORMP  perform  a polynomial  interpolation  using  HERMITE'S 
interpolation  formula,  INTGRLN  and  INTGRLP  integrate  the  area  under  the  curve 
between  given  limits  of  A and  B. 

Programs  HERMITN  and  HERMITP  construct  the  polynomial  using  the  same 
HERMITE'S  interpolation  formula,  and  CORM4  performs  polynomial  interpolation 
using  Lagrange's  inter}  iation  formula.  The  difference  between  the  programs 
of  GROUP  I and  GROUP  II  is  that: 

GR0UPI:  Constructs  a polynomial  going  through  the  points  having 

positive  or  negative  value. 

GROUP  II:  If  the  evaluation  of  the  constructed  polynomial  becomes  negative 

(that  is  unreal  for  wind  speed  data)  then  the  programs  of 
GROUP  II  correct  the  negative  values  and  such  that  the 
polynomial  evaluation  is  always  positive. 
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To  accomplish  the  tasks  outlined,  the  following  mathematical  steps  have 
been  taken: 


A.  Program  CORMN  and  CORMP  perform  a polynomial  interpolation  by  using 
HERMITE'S  interpolation  formula.  Suppose  that  values  of  f(x)  and  f ' (x)  are 
known  for  x1,...xm.  A polynomial  H(x)  can  be  determined  by  assuming  that  it 
is  expressible  in  the  form: 


m m 

H(x)  = E h (x)  f(x.  ) + E H.  (x)  f'(x.) 
k=l  K K k=l  K * 


where  h^(x)  and  h^(x)  (i=l,2, . . . ,m)  are  polynomials  of  maximum  degree  2m-l: 

h^x)  = [l-2L[(x.)(x-x.)][L.(x)]2 


^(x)  = (x-xA)  [L±  Cx)  ] 


where 


Li 


(x-Xj) (x-xi_1)  (x-xi+1) . 

(x.-x1y:...(x.-xi_1)  (x'i-x^p 


(x-xn) 


^vv 


For  two  points  we  have  a polynomial  of  degree  three: 


H(x)  = [f(x1)V1(x)+f» (x1)W1(x)]L12(x)+[f(x2)V2(x)+f»(x2)W2(x)]L22(x) 


where 


x-x„ 


L.  = — V = 1 

1 X,-X„  1 


'1  2 


(x-x1)u)"(x1) 
u ' (Xj) 


= x-Xj 


x-x. 


V.  (x-x-)w"(x-) 

I - 1 V _ i - f 

L2  “ x^-x,  V2  u* (x7) 


'2  1 


W2  = X-X2 
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(1) 


(2) 


* 


-■  » i-  ■ ■ ■ - -- 


and: 


U(x)  = (x-Xj) (x-x2) 


U)'(x)  = (X-Xj)  + (x-x2) 

“'<xi)  • xrx2 

w'(x2)  = *2~xl 
u"  = 2 


B.  Programs  INTGRLN  and  INRGRLP  integrate  the  function  (2)  for  any  given 
limits  of  A and  B under  the  curve  H(x) . The  following  steps  have  been  taken: 

ANS  * /fCxpVjWL^Mdx  ♦ /f(x2)V2(x)L22(x)dx  ♦ ft ' (x^WjCxJL^dx 

♦ ff' (x2)W2(x)L22(x)dx 


where 


Y1  = f(Xl) 


Y2  = f(x2) 


SL1  = f’CXj) 


SL2  = f'(x2) 


PI  = Xj-X2 


P2  = x2-Xj 


are  given. 

By  integrating  each  term  of  the  equation  (3)  we  get 

/f(x1)V1(x)L12(x)dx  = fCx1l/V1(x)L12(x)dx  = f(Xj) 


(x-x. ) 2 x-x_  2 
/(I  - -p--  ) (-pr— ) dx 
*1  *1 

£(Xi)  2 

dx  + *— ./*(x-x,)  dx 

V 


2f(xJ 


■^-/(x-Xj)  (x-x2)‘ 
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y y 

/f(x2)V2(x)L2  (x)dx  = f(x2)/V2(x)L2  (x)dx  = f(x2)/(l p ■) 


x-x.  2 m j 

(-p— 3 dx  3 — /(x-x2)(x-x1)  dx  ♦ — 2 /(x-x/  dx 


2f(x^) 

P. 


f(x,) 
P„ 


(3b) 


(X-X  ) f • (x  j 

/f'(x1)W1(x)L12(x)dx  = f'M/U-Xj)  y — dx  = ji-^tx-Xj)  (x-x2)2  dx  (3c) 

P P1 


(X_X  j (X  j 

/f'(x2)W2(x)L12(x)dx  = f ' (x2)/(x-x2)  5 — dx  = ^-Jlx-x^  (x-Xj)2dx  (3d) 

P2  P2 


Since: 


F = /(x-a)  (x-b)2dx  = j(x-b)3(x-a)  - y/(x-b)3  d(x-b)dx  = i(x-b)3 


(x-a)  - = F(a,b,x) 


and 


G = /(x-a)2  dx  = S&L-  = G(a,x) 


For  the  given  limits  a and  b of  the  integration  Equation  (3)  becomes: 
ANS  = (C1+E1)*(F1-F2)+(C2+E2)*(F3-F4)+D1*(G1-G2)+D2*(G3-G4) 

where 

Cl  = -2.0*Y1/(P1**3) 

C2  = -2.0*Y2/(P2**3) 

D1  = 1.0/(P1*P1) 

D2  = 1.0/ (P2*P2) 

, \ 
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= SL1*D1 
= SL2*D2 
= F(X1,X2,B) 

= F(X1,X2,A) 

= G(X2,B) 

= G(X2,A) 

= F(X2,X1,B) 

= F(X2,X1,A) 

= G(X1,B) 

= G(X1,A) 


- f(xx) 

- f(x2) 

= f'Cx^ 
= f'(x2) 


xrx2 


x2'xl 


C.  Program  C0RM4  performs  the  Lagrangian  interpolation  (given  N points  to 
fit  exactly  by  a polynomial  of  degree  n-1)  using  the  formula: 


£ a y./(x-x.) 

j=1  l i 

n 

r a./(x-x.) 


where  the  given  data  are  (x^y^,  i=l,2,...,n 


aj  " (x^-Xjj  (Xj-x2)'.V.  (xj-Xj  j)  (xj+1)".V/txj-xn) 


I nveteed  Depth  to  Total  Depth 

Initiator:  Rene  V.  Cormier 

Problem  No:  4026-7  Project  No:  8624 


The  purpose  of  this  problem  under  the  title,  "Inverted  Depth  to  Total 
Depth",  is  to  study  the  meteorology  of  integrated  boundary  layer  winds  in 
support  of  Air  Force  paradrop  operations.  Specifically,  the  effect  of  atmospheric 
stability  on  integrated  boundary  layer  winds. 

Conventional  methods  of  determining  atmospheric  stability  are  inadequate 
to  define  the  stability  of  the  atmospheric  layers  under  integration.  This 
program  develops  a new  stability  parameter.  It  considers  the  temperature 
structure  throughout  the  depth  of  the  layer  rather  than  just  top  and  bottom 
temperatures.  It  finds  the  ratio  of  the  depth  within  which  the  temperature  is 
constant  or  increases  with  height  to  the  total  depth  under  integration. 

The  following  steps  have  been  taken: 

Polynomial  interpolation  is  performed  on  the  input  vertical  temperature 
profile  using  Hermite's  interpolation  formula  in  order  to  perform  the 
integration. 

The  coefficients  of  a defining  cubic  polynomial  function  are  evaluated. 

The  maxima  and  minima  of  the  function  are  then  found  by  considering  the 
first  and  second  derivatives  of  the  function. 

To  accomplish  the  tasks  outlined,  the  following  mathematical  steps  have 
been  taken: 

The  polynomial  interpolation  and  integration  is  done  using  the  method  and 
procedures  described  in  problem  number  4026-6,  project  8624  above. 

Evaluate  the  coefficients  of  the  cubic  polynomial  H(x)  = ax3+bx2+cx+d. 
y!+y'  2(y  -y  ) 

a = — 2 — + ! — i- 

<vv2  (vv3 
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From  the  two  point  Hermite's  interpolation  formula 


Expand -con  o £ GRAPPAC2 

Initiator:  Mr.  R.  Gosselin 

Problem  No:  4040-1  Project  No:  0002 

CALCOMP  plotter  simulation  capability  on  CDC's  Interactive  Graphics 
System  was  made  available  to  all  computer  users  through  expansion  of  the  GRAPPAC 
system  to  include  all  CALCOMP  entry  points.  An  instruction  booklet,  skeleton 
deck,  and  associated  permanent  files  were  prepared.  As  a result,  CALCOMP  jobs 
can  now  be  debugged  quickly  and  easily  without  waiting  for  hard  copy  test  plots 
to  be  prepared. 

In  addition,  CALCOMP  type  programs  can  now  be  made  to  accept  input  from  * 
the  graphics  console  through  the  addition  of  a very  simple  subroutine  call. 
Selective  erasure  is  available. 


GRAPPAC:  A package  of  Fortran  Subroutines  for  use  with  the  6000  Series  274 

Interactive  Graphic  System  of  the  Control  Data  Corporation,  Frona  B.  Vicksell, 
Space  Data  Analysis  Lab.,  Boston  College,  Chestnut  Hill,  Mass.  02167,  AFCRL-72- 
0698. 


A naly&iA  orf  Ob&zxved  Election  VejuiXy  Pio^ilu 

Initiator:  Mr.  R.  Allen 

Problem  No:  4041-1  Project  No:  8666 


Support  was  provided  for  study  of  a set  of  observed  electron  density 
profiles.  The  profiles  were  obtained  from  Lincoln  Laboratory,  Bedford,  Ma. 
Programs  were  written  first  to  present  the  data  in  graphical  form  for  visual 
screening,  then  to  perform  certain  calculations,  producing  a new  data  base 
for  use  in  future  statistical  studies.  The  calculations  for  each  profile 
involved  integrating  the  available  densities,  with  extrapolation  if  needed, 
over  the  altitude  range  200  to  1000  km. , to  obtain  total  content  and  slab 
thickness.  In  addition.  Chapman  functions  and  parabolas  were  fitted  to  the 
peak  density  region  to  obtain  estimates  of  the  density  scale  height  above 
the  peak,  below  the  peak,  and  on  both  sides  simultaneously. 


Modeling  Election  Ven&iXy  Pno^itu 

Initiator:  Mr.  R.  Allen 

Problem  No:  4041-2  Project  No:  8666 


In  a study  related  to  account  no.  4041-2,  statistical  and  graphical  com- 
parisons were  made  between  three  diffeqent  sets  of  observed  electron  density 
profiles  and  a model  provided  by  the  Air  Weather  Service  and  modified  by  the 
Initiator.  The  model  has  a variable  scale  height  Chapman  function  on  the 
topside  of  the  F2-region,  a fixed  scale  height  for  a distance  of  one  scale 
height  for  a distance  of  one  scale  height  below  the  F2  peak,  then  linear  inter- 
polation on  the  logarithms  of  the  densities  down  to  the  E-region  peak.  The 
bottom  side  of  the  E region  is  another  Chapman  function. 

Inputs  to  the  model  include  the  peak  F2  density  and  its  height,  latitude, 
longitude,  date,  time,  and  10.7  cm.  solar  flux  or  sunspot  number. 


Theoretical  Prediction*  the  Electron  Vi&tribution 
Function  in  the  Ionosphere 


Initiator:  Dr.  J.  Jasperse 

Problem  Nos:  4585,  4688,  4799  Project  No:  6688 

These  problems  involved  making  theoretical  predictions  of  the  electron 
distribution  function  in  the  ionosphere.  Over  a period  of  time,  the  model 
considered  has  been  made  more  and  more  extensive  until  the  following  integral- 
differential  equation  was  solved  for  H(z,E)  in  the  latest  computer  program. 


2 in 

> w*>e>  e’/2  (=-=-)  yi..w 


i . +m  'mnj 
nj  e J 


'm.  +m  ' 'mis 
is  e 


♦ K I°(z 

3 0 


,E)  I H(z,E)  =J  2dEJ  ^ G^Cz.E1; 


d£l  ^ QlCz.El)  j*2  ^ (E1) 1/2  Y^Cz.E1)  HCz.E1) 


2 -E^TW* 


(E1)1'2  Y^Cz.E1)  6 (E*-E  , ) Hfz.E1) 


+ I dE 


2 


(E‘*Etk>  W2-E‘*Etk>  H<2*El*Etk> 


dE‘  e’‘P  1-  icTTHTl  (E‘*Etk)‘/!  l'tkt2'E‘*Etk>  H<2-E‘> 


MJOt) 


^ t~^TntO]  (£l)l/2  Ytk(z*El)  0(El-|W  H^.E‘-Etk) 


WJ 


Z^j  <v:%)  W‘-E)  E ',T«w 


NS 

*5  <=£%>  ^ W*> 


+ 3 ° C*» E)  + j (z»E)]  1 H(z,E) 


Where: 


W**B>  = % <*>  (CE)l/2  <^(B) 


Ymis^z>  " K2  Nis<2> 


YrJlCz,E)  * N.s(z)  (CE) 1/2  Q^CE) 


Ytk(2*E)  = Nnt(z)  (CE)  Qf-vCE) 


lJ(*.E)  = / dE1(E1)l/2  H(z,E!) 

•'P 


(El) 3/2  H(z,E‘) 


I 

l 


A 

1^(2, E)  = E3/z  / 2 dE1  H(z, 

J p 


E1) 


MMjJfJ) 

vz'e)  ■ y,  vz)  Vi 

m=l 


ljm(E*Epij«>  I(J’E*Epij») 


I(z,E)  = I (b,E)  exp  1 -sec  x 


Vi(E)  / dz‘  Vz‘> 


Gi(z-E,)  ■ Ks  %® 


£/> 

K i t +ti;jk 


dE11EuD.jk(Ell,E1)  H(z,Eu) 


- E1  H(.#E»)  E e(B*.Eijk)  Q..k(E3)] 


etx-x1)  = 


0 x < x1 

1 x > x1 


n re11  pi)  _ qo  ijk  ijk  s 

iiklfc  J 7, r <— VT  ' — ) 


E. 


'ijk' 


2(E l/z+E.  ..)*  ''ei/2+E. 

ijk''  c ijk 


■ li 


. ^"nkj  °ijk  n . tEi/2^i1t.Btfi  v 


-li 


-) 


■] 


For  E...  < E11  < E and  0 < E1  < Eu-E.., 
ijk  — — 2 — — ijk 


Dijk(Ell»El)  = 0 for  all  other  values  of  E1  and  E11 
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< ..  ~ 


<1EU  Dijk(E1’Ell) 


For  t»l,  k-1,7  MK(1)=7 


Qtk(E) 


0 for  E.  < E < E., 

1 — tk 


fqo  Atk.  ,Etk 


Etk,  V\k 


) c-r->  °tk  ti  - (-£)  ) 


For  E..<  E < E 
tk  — — 2 


For  t=*2  MK(2)=23 

For  t=2,  K»l,  8 Q2K(E)  is  linearly  interpolated  from  a set  of  Q K(e)  versus 
E values 


For  t=2,  K=9-23 


0 for  E < E < E , 
l — tk 


( Etk 


For  E..<  E < E 

lK  — — 2 


For  t=3  MK(3)=14 


For  t=3  K=l-4  Qjk(E)  is  linearly  interpolated  from  a set  of  QJk(E)  versus 


E values 


For  t*2  K=*5-14 


I 


0 For  Ej  E < E^ 


*tk 


<\n  Kv  E,v  °tk 


C-jp)  11 


'tk 


Ak  “tkMtk 

(— ) ] For  Etk  < E < E^ 


For  t=4  MK(4)=2 


For  t=4,  K=l,2  Q , (E)  is  linearly  interpolated  from  a set  of  Q .(E)  versus 

4 K 4 K 

E values. 


For  t=5,55  MK(t)=l 


For  t=5,K=l  Q (E)  is  linearly  interpolated  from  a set  of  Q (E)  versus 
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E values 


For  t=6-35  K=1 


Etl  = BjCt-4)  (t-3)  - (t-6)  (t-5) ] 


and 


| ; 


MntU)  • V*) 


(2t-71)  exp[-  gT^-] 

; (2U-71)  expt-  r^y] 


• = q2 


K 


dE  E Q^fE)  H(z,E) 


(K  +K  ) Nn  + CK  +K  ) N + K y N _1 

211  214  111  233  2 34  n3  244  ' n3 


[q  + K 
1 211 


V vjc,/*  / 


dEE  Q^fE)  H(z,E) 


+ K N + K N + K Y N I"1 
124  n2  l33  n3  144  n 3 ( 


q + K N N.  + K N N.  ] 

3 2 3 3 n3  12  133  n3  liJ 


Cl/z  / 2 dE  E 0(E)  H(z/E)  + K N + K y N 

J ”3  324  n2  344  H 3 

E 


N+K  N +K  YNlN. 

214  ni  2 34  n3  244  T Il3J  12 


Nn  + K Y ] N.  + [K  N + K YN  In. 

124  I>2  144  n3J  11  L 324  n2  344  1 Il3J  l3j 


cl/z  / 2 dE  E Qy^E)  H(z,E) 
•'e 


i 


is  approximated  as 


For  £ = 1 = H(z,EP  ^ ) * 2. 248. 10'2  0 ( (EMID^+1))  3/2  - (EMID^V^2) 

For  £ = 2 = H(z,EP(v))  * 6. 426727273. 10'1 6 ( (EMID^ +1^ ) * * - (EMID^V  *) 

For  £ = 3 = H(z,EP^)  * 3. 033666666710"16  ( (EMID^+1))  * 9 - (EMID(^)’9) 

For  £ = 4 = H(z,EP(vj)  * 4. 39163636364. 10" 16  ((EMID(;i  + 1) ) 1 * 1 - (EMID^)1’1) 


where  EP^  is  the  EP  value  such  that 


EMID^  lEP(v)  iEMID^+1‘) 


K221  = K211S 


K214  = K214S  Cy^y)  4 ♦ ^214SS 


1000  "8 

K233  = K233S  (y^) 


1000  ' 4 

K133  = K133S  Cy-yyy) 


Tn(z) 

K144  = K144S  (ySoo“) 


K124S  (sr^^V)  for  T < 750 
Tn(z)'  n - 


Tn(z) 

K124S  (-750-)  for  Tn  > 750 
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For  these  calculations  we  have  the  following  input  data 

K211S,  K214S,  K214SS,  K233S,  k234,  K244,  K124S,  K133S,  K144S, 

K324,  K344,  GAMMA,  M . K , K , K , K , K , X,  C,  a set  of  EP,.. 

e l 2 3 H 5 llj 

values,  a set  of  Z values  (Z(i)),  NJ,  Mnj,  qQ,  NK(J),  Eijk,  A.jk, 

Pijk*  ♦ijk*  Bijk*  Nijk’  NS«  Mis*  VZ)’  Nn,*(Z)’  NL 

(NL  = NS  is  implicitly  assumed),  MW(J),  E . . , E . , A . , <p  . , B . , 

pxjiu  1 lv  IK  l 

Eik’  Ajk*  *,k’  B;k'  N2k’  E,k-  Aik-  ♦jk’  B,k’  Njk’  E,,'  E»2-  E,,>  qi‘  Eci* 

V V EC2-  B2 


In  addition  the  program  reads  the  following  bar  graphs: 
Bar  graphs 

NI.  values  of  I.  versus  E 

D D 

NQPAJj  values  of  versus  E 
NQPIJKj^  values  of  Qpjjm  versus  E 


And  finally  the  program  reads  the  following  graphs  which  are  linearly 
interpolated 


NTN  values  of  T versus  z 
n 

NT  values  of  T.  versus  z 

i 

N(»1t  values  of  P . versus  z 
J mi 


Given  the  preceding  data  the  differential-integral  equation  for  H(z,E) 
is  approximated  by  a difference  equation  for  a particular  value  of  z. 
To  do  this,  the  equation  is  evaluated  at  the  points  EP^. 
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The  integral  part  of  the  equation  is  usually  approximated  by  assuming 
H(z,E)  is  a bar  graph  with  values  H(z,Ey)  covering  intervals  [EMID^, 
EMID(-L+1j].  Exceptions  to  this  rule  are  for/dEEQre  H(z,E]  and 

/dE‘S  GjU.E1)  which  are  handled  as  previously  described  (for  fd\l  L Q 

H(z,E)).  Another  exception  is  that  above  a read-in  value  for  E the  terms 


- -f  * dE' 


Mjcm 


(El) 1/2  Ytk(z,E‘)  0(El-Etk)  HCz.E1) 


/ a 

+ / dE1 
E 


(E‘*Etk)‘/2  Ytk(l,E‘.Etk)  H(Z,E>*Etk) 


“f  I-  irrm-l  y.Jj.e'.e  .)  hu.e1) 


KTn(z) 


E 55  NIK  ft]  i 

{ dE‘2  S 6XP  ["  KTTi7>  (E‘>‘/2  Y^Ci.E1)  6CE‘-Etk)  H(Z,E'-B  ) 

f*=  7 ^^7  1 n 


55  MKft) 


2l  11  Etk exp  [-  rrriy]  t<VEtk>l/2  W-VW  H(2*V 


t»l  k=] 


(E+Etk),/2  Ytk(z.E+Etk)  H(z,E)] 


Respectively  where  we  have  used  the  approximation 


- fCE-Etk)  _ , f (E 


d HE) 

W • 


The  derivative  part  of  the  equation  is  found  at  the  values  EP^  in  terms 
of  central  differences  generally  a seven  point  rule  is  used  (a  smaller  odd 
point  rule  can  be  read  in  as  this  input  parameter)  except  near  EP^  and 

EP (NE) * Thus>  Ep(s)  and  EP(NE  2)  use  Eive  Point  ru1®5*  EP(2)  and  EP(NE-i) 
use  three  point  rules  and 

■ 3 H(z,E)  __  H(z,EP(N£)  - H(z,EP(Nfii)) 


(NE  " (NE-,) 


E = EP. 


The  equation  at  EPf  . is  replaced  by  the  equation 


= dE(E)l/2  H(z,E)/Ne(z) 


1^— 2*-E--  is  never  approximated. 


E * EP. 


Using  the  logic  specified  above,  the  problem  was  reduced  to  the  solution 
of  NE  non-linear  equations,  one  for  each  EP^  value  for  the  values  of 


I 


H(z,EP^yj)  for  a specified  z value.  This  set  of  equations  was  solved  for 
H(z,EP^)  using  a modified  Newton-Raphson  method.  Thus,  if  AH(z,EPri^) 


(v)' 


is  the  correction  calculated  by  applying  the  Newton-Raphson  method  to  the 

NE  equations.  Then  the  corrected  H(z,EP^)  = the  old  H(z,EP^)  + 

2~k  AH(z,EPy)  where  L is  the  smallest  integer  which  reduces  the  value  of 

NE 

ET<J>T.  Where  ET$T  = SCA(*EREEQ(1))2  + ^ (EREEQ(i))2.  Where  EREEQ(j)  is 

the  amount  by  which  the  equation  is  not  satisfied  when  the  guessed 
values  of  H(z,EP^)  are  tried.  The  program  continues  to  correct  H(z,EP^vj) 
until  a "solution"  is  found. 


In  addition  to  the  above  solution  the  program  calculates  and  plots  the 
following  quantities  when  asked  by  the  input  data. 


I 

!!- 

r * 

M 

( 


Oxygen  Excitation  Cross-section  (versus  Energy) 


Q K(E)  + Q (E)  + Q (E)  + Q (E) 

1 41  42  51 


Oxygen  Ionization  Cross-section  (versus  Energy) 


QiiK<E> 


N Rotational  Cross-section  (versus  Energy) 
2 


Qtl(E) 
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Excitation  Cross-section  (versus  Energy) 

■ 

N Ionization  Cross-section  (versus  Energy) 
2 

0 Rotational  Cross-section  (versus  Energy) 
2 


0 Excitation  Cross-section  (versus  Energy) 
2 


0 Ionization  Cross-section  (versus  Energy) 
2 


I 


(E) 


After  each  iteration  of  H(z,EP^)  the  program  calculates  and  points 


the  following  quantities: 


fE 

2 dE  F(z,E)) 


where  F(z,E)  = y/fi  H(z,E) 


f2  E*F(z,E) 


T = 7736.666667  (I  2E*F(z,E)  dE)/N  (z) 
e E e 

l 


M2>  * ^ Ni, 


Integral  of  G1  = J^2  dE  \^Gj(z,E) 

Ei 


and  each  integral  of  Gj 


= J 2 dE  G?  (z,E)  for  J = 1,2,3 


When  the  program  finds  a solution  for  H(z,EP^),  the  program  also 
calculates  and  prints  (plots) 


F1 (Z»E)  = 


' [1  - exp(.  r^T)J  Etk  YtkC*.E)  e(E-Etk) 


NJ  2 m MS  2 m 

* E E(iw5:)  ymi  * E "7? 

J nj  e J VE  ij  e 


♦ — V 1 (z)  K El/z 
y/E  66  5 


Eijk  VZ)  0CE-Eijk)  Qijk(E) 


C dE‘  Sgj(i,e1) 
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Check  for  G ss  Vz)  /*  dE1  E1  Qi^CE1)  H(z,E») 


ijk 


conductivity 


6.  a T K f2  dE  E3/z  L (H(z,E))/f 

E oz  j 


"13  T 4 


Wz>«  * rl»is(l)) 

S*1 


approximate  conductivity  */•  (z) 

<Ss  'approx. 


c / 

t *D 


■ (S^(e-,,f^’/!Sy^w) 


Average  Ionization  Collision  Frequency  = ^(z) 


(to-) 


^ NK^  £ 

/Ai  (z)  T!  / 2 dE1  (E1) 1/2  Qijk(El)  FCz.E1) 


ljTJzT'  — — -F 

j=l  k*l  Eij 


Alternate  calculation  for  Electron  Concentration 


SS 


j 


Ni 


V f iL(Z)  V'J>) 


NL  P NJ 

t4  Z^fiL(z)  “rL(z)  { 1 d£l  J^G.Cz.E1)  ♦ V*(z)]  ♦ V. 
L=1  1 j=l 


U) 


arL(z)  for  L = 1,2,3  and  4 


■ rr^r  f‘ dE  E 


ev  ' e 


Approximate  Average  Ion  Collision  Frequency  = 


K ^ NK(J)  p 

Vz)  = N“hrl-/Nnj(2)  S J 1 dEl(£l)l/2  Qijk(El)  pl(z*El) 

6 j=l  K=1  Eij 


Approximation  for  the  Electron  Concentration  Using  Vp(z) 


mP  _ r 1 

e - [f  T y . ■ - , ^ 

L=1  fiL(z)  arL(z) 


•] 


4 > . fiL(zi  arL(z) 


A 

f dE1 


G.(z,El)  + Vp(z)]  ♦ Vp(z) 


Maxwellian  Approximation  to  F = F^fz.E) 


4lT  Ne(z)  ( 2 } [wk  t czy]  exp(_  rr(z)} 

i e'  j e*-  J 
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Early  in  the  program  the  program  can  read  in  guesses  for  N (z),  N.  , 

e lx 

Ni2»  Ni}.  and  Te(z)  and  use  (FM(z,EP(y))  + F1 (z,EP(v))/(EPIv))l/2  as 
the  initial  guess  for  Hfz.EP^)  to  start  the  iteration  procedure. 


Finally,  for  a given  z value  the  program  calculates 


Ex  = 


(CE) 


1/2 


55  MK(t)  E 

EE11' exp(-  irrW>  Etk  *tklI-E>  e<E-Etk) 

t=l  k=l  1 " 


NJ  2b  NS  2 m 

Ee  W*-»  * (s-s-)  *ii,w 


3*1 


nj  e 


s=l 


(E)1 


13  e 


NJ  NK(J) 

* * K,  E'/2  E E Eij  V2)  e(E-Eijk)  «ijktEij 

1 J j=l  k=l 


This  program  will  be  submitted  to  the  SUYA  computer  program  library 
upon  completion  of  problem  4799.  The  Initiator  intends  to  further  refine 
the  equations  for  H(z,E).  In  any  case,  a separate  program  has  been  written 
which  uses  the  calculated  H(z,E)  functions  to  calculate 


SF(z,E)  = 4.72* 10s  E H(z,E) 


4 

4 «. 


SF1 (z,E)  x 4.72- 106  E H‘(z,E) 


qex(z)  = S.93-107  N (z)  jf  6X2  dE  Qq(E)  E H(z.E) 

exi 
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where  Nq(z),  fio(z)>  T(z),  NgCz),  E , EeXl  are  input  parameters  and  Qq(E) 
is  an  input  function. 

The  results  of  the  program  seem  to  agree  with  experiments  for  low  values 
of  altitude,  up  to  about  200  kilometers.  In  fact,  they  predict  fine  structure 
which  the  instruments  which  have  been  used  to  measure  electrons  in  the 
ionosphere  were  too  coarse  to  find.  Above  about  200  kilometers  the  model 
needs  improvement.  The  initiator  hopes  to  refine  the  model  even  further 
and  explain  the  electron  distribution  above  200  kilometers. 
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Study  Sous ice  Location  by  Compute si  Simutation 

Initiator:  Dr.  K.  Toman 

Problem  No:  4620  Project  No:  5631 


Locating  sources  of  radiation,  finding  radar  targets,  and  determining 
one's  own  position  are  requirements  in  radio  surveillance,  detection,  and 
navigation  respectively.  There  also  exists  in  the  environmental  sciences  the 
need  for  determining,  for  example,  the  location  of  the  origin  of  natural  or 
man-made  atmospheric  or  seismological  disturbances. 

In  the  following  treatise  we  have  restricted  ourselves  to  study,  on  a 
computer,  the  problem  of  locating  sources  from  time-of-arrival  measurements 
made  at  four  locations.  We  developed  and  tested  an  error  analysis  to  determine 
how  errors  in  time  differences  reveal  themselves  as  errors  in  source  location 
and  signal  velocity.  Furthermore,  we  ascertained  how  these  errors  change  with 
range,  azimuth,  signal  velocity,  and  network  configuration.  The  results, 
subject  to  certain  assumptions,  can  be  used  for  existing  or  planned  time-of- 
arrival  measurement  networks.  These  assumptions  are:  (a)  source  and  synchronized 

receiving  stations  are  located  at  the  surface  of  a spherical  earth  along  which 
the  signal  propagates,  (b)  the  velocity  of  the  signal  is  uniform  along  the 
great-circle  paths  from  source  to  receiving  stations,  (c)  the  signal  is  not 
bandwidth  limited  nor  is  it  dispersed  in  frequency,  and  (d)  individual  errors 
are  normally  distributed  and  uncorrelated*.  Assumption  (a)  can  be  broadened 
to  include,  in  an  idealized  form,  the  case  of  remote  sensing  of  atmospheric 
waves  by  ionospheric  sounding  techniques.  Assumption  (b)  can  be  broadened  to 
include  rectilinear  propagation  and  reflection  of  a radio  signal  from  a mirror- 
like ionosphere  of  known  height. 

The  locations  of  the  source  S and  of  the  receiving  stations  A,  B,  C,  and 
D on  the  surface  of  the  earth  of  radius  R = 6378  km  are  expressed  in  geographical 
coordinates  0 (longitude)  and  <p  (latitude) . Using'  differences  At  of  the  arrival 
times  for  the  propagating  signal  at  the  synchronized  stations,  the  location  of 
the  source  (0g,  <f>g)  and  the  signal  velocity  (v)  can  be  determined  from  the 
following  equations: 


* Cooper,  D.C.  and  Laite,  P.J.  (1969)  Statistical  Analysis  of  Position  Fixing 
in  Three  Dimensions,  Proc.  IEE  116  (No.  9) :1505-1508. 
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A~B  — 

= cos'1  [sin<J>A  sin<J>g  + cos<pA  cos<f>g  cos(0s-0A)] 


cos'1  [sin<f>  sin0  + co s<p  cos<p  cos(0  -0  )] 


Y~~  s cos'1  £sin4>c  sin<f>g  + cos<l>c  cos<|>g  cos(0g-0c)] 

- cos'1  [sin<j>B  sin<J>g  + cos<j>B  cos<f>g  cos(0g-0g)]  . 


k — = cos'1  [sin<J>D  sin<J>s  + cos<J>D  cos<J>g  cos(0s-0D)] 

- cos'1  [sinc^g  sin<J>g  + cos<J>B  cos<t>g  cos(0g-0g)]  . (1) 

In  the  method  used  for  obtaining  the  solutions  0g,  <f>g,  and  v from  the  time  differ- 
ences AtA_B>  Atc_B,  and  AtQ_B  of  known  arrival  times  tA,  tg,  tc,  and  tD  of  the 
signal  relative  to  a common  time  reference,  two  of  the  equations  in  Eq.(l)  were 
combined,  a speed  v^  was  assumed,  and  0g^,  <J>g^  were  determined.  These  "initial" 
source  coordinates  were  entered  into  the  third  equation  to  determine  a speed 

v which  was  compared  with  v. . If  v-  < v then  v.  was  increased  and  vice  versa, 
c 1 X c X 

This  iteration,  which  employed  the  False  Position  Method^,  was  continued  until 
the  difference  between  v^  and  v£  was  negligibly  small  (for  example, 

| vi  vc  1 - 1 <10  6).  The  initial  assumption  for  v^  was  thereby  only  limited  by 
considerations  of  minimizing  computer  search  time.  The  iteration  procedure  was 
used  in  the  neighborhood  of  v=18000  m/min.  For  a given  source  location,  arrival 
times  at  the  four  stations  and  their  time  differences  were  computed  using  this 
speed.  For  performing  the  inverse  task,  these  time  differences  were  used  to 
redetermine  the  source  location  and  velocity  prior  to  assuming  that  the  former 
are  subject  to  errors.  In  order  to  determine  the  effect  of  these  time- 
difference  errors  on  errors  in  source  location  and  signal  velocity  an  error 
analysis  was  developed.  Errors  were  imposed  upon  time  differences  At  by 
assuming  that  they  represent  standard  deviations  of  a normal  distribution.  The 
resulting  errors  in  source  location  and  signal  velocity  were  viewed  in  terms 
of  probability  bounds. 

2 

Scarborough,  J.B.  (1966)  Numerical  Mathematical  Analysis,  John  Hopkins  Press, 
Baltimore,  6th  ed.,  p.  197. 
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In  this  error  analysis  it  was  assumed  that  the  errors  of  the  independent 
variables  AtAB,  Atc_fi,  and  Atp  B of  Eq.(l)  are  small  and  uncorrelated.  Based 
on  these  assumptions  the  formula  for  the  propagation  of  the  mean-square  error^ 
could  be  used  to  estimate  the  mean-square  errors  OQg,  o|s,  and  of  the 
dependent  variables  0g,  <j>g,  and  v as  follows: 

_2  /39S  \2  /39s  \2  /36s  \2 

es  ■ WTe  X-bJ  * (^0N-b)  • 

2 

= (Replace  0g  in  the  above  with  <f>g)  , 

S 


(2) 


= (Replace  0g  in  the  above  with  v)  , 


; 


where  , ...  etc.,  represent  the  standard  deviations  from  the  mean  assumed 

A-B 

for  the  time -of- arrival  differences  between  selected  station  pairs.  As  seen 
from  Eq.(l),  explicit  expressions  for  0g,  <|>g,  and  v were  not  available.  The 
evaluation  of  their  partial  derivatives  30g/3AtA  B,...  etc.  [Eq.(2)]  required 
the  use  of  a special  property  of  the  Jacobian  matrix  and  its  inverse2*. 


i 


P4tA-B 

9AtA-B 

“VbI 

-1 

r96s 

"s 

"s  1 

"s 

S+S 

3v 

9AtA-B 

9AtC-B 

9AtC-B 

3itC-B 

34tC-B 

3^S 

C/D 

3*s 

3+S 

3v 

3itA-B 

9AtC-B 

9AtD-B 

34tD-B 

34tD-B 

3md-b 

3v 

3v 

3v 

"s 

3v 

5^ 

34tC-B 

^D-B 

(3) 


whereby  the  partial  derivatives  3AtA_B/30g, . . . etc.,  were  obtained  from  Eq.(l). 
Using  this  approach,  the  mean-square  errors  in  0 , 4>  and  v were  computed 

S b 


Deming,  W.E.  (1964),  Statistical  Adjustment  of  Data,  Dover  Publications, 
New  York,  p.  39. 


Irving,  J.,  and  Mullineaux,  N.  (1959),  Mathematics  in  Physics  and  Engineering, 
Academic  Press,  New  York,  p.  800. 
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by  means  of  Eq. (2)  from  the  errors  assumed  for  the  time  differences.  The 
results  of  these  computations  were  used  as  error  bounds  (error  boxes)  centered 
on  true  source  locations.  These  error  bounds  were  tested  by  generating  a 
sequence  of  independent  errors  in  time-of-arrival  differences,  using  distribu- 
tions for  means  and  variances  that  were  derived  from  available  random  number 
routines.  Although  these  bounds  were  found  to  contain  the  appropriate  number 
of  "false"  source  locations,  it  was  noted  that  their  distribution  seemed  to 
fill  these  error  boxes  only  partially,  due  to  an  elongated  scatter  of  points 
in  a direction  toward  the  network  with  relatively  little  scatter  in  azimuth. 

It  was  also  noted  that  the  inherent  nonuniformity  of  the  geographic  coordinate 
system  would  give  the  appearance  of  large  errors  in  terms  of  oqs,  when  the 
source  is  near  a geographic  pole,  and  the  appearance  of  small  errors  when  it 
is  near  the  equator.  For  these  reasons  it  was  decided  to  describe  the  "probable 
location  of  the  source"  not  only  in  terms  of  longitude  and  latitude,  but  also 
with  respect  to  distance  and  azimuth.  This  made  it  necessary  to  choose  a 
reference  point  relative  to  which  distance  and  azimuth  are  defined.  Although 
any  point  on  the  surface  of  the  earth  could  be  chosen,  one  of  the  stations,  B, 
of  a network  was  usually  designated  the  reference  point. 

I 

Equations  were  obtained  that  relate  the  coordinates  of  a source  in  latitude 
(4>g)  and  longitude  (0g)  to  range  (Dist)  and  azimuth  (Az)  for  a reference 
station  B C<I>B,  9g) . 

<J>S  = sin-1  | sin  4>B  * cos(Dist/R)  + cos  (J>B  • sin(Dist/R)  • cos  Az  j 

0 

9g  = 6g  - sin  ^sinfDist/R)  • sin  Az/cos  <J>S  j.  . 

From  these  equations  the  partial  derivatives  34>„/3Az,  3<J> „/3Dist,  30  /3Az, 

v)  b O 

30^/3  Dist  were  determined  and  related  to  the  partials  of  the  time  differences 
At  with  respect  to  azimuth  and  distance  by 
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where  3At/30<.  and  3At/3<f>g,  as  well  as  3At/3v  were  obtained  from  Eq.(l).  The 
partial  derivatives  of  azimuth,  distance,  and  velocity  with  respect  to  the  time 
differences  that  were  to  provide  distance  and  azimuth  error  boxes  analogous  to 
Eq.(2)  were  obtained  from  the  partial  derivatives  of  the  time  differences  with 
respect  to  azimuth  and  distance  [obtained  from  Eqs.(4)  and  (5)]  as  well  as 
velocity  [obtained  from  Eq.(l)]  by  utilizing,  as  before,  the  Jacobian  inverse 
property. 
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The  above  equations  were  used  to  evaluate  the  errors  of  source  location 
(<J>S>  9g,  or  Dist,  Az)  and  signal  velocity  v from  given  errors  in  the  time 
differences  and  to  appraise  the  effects  of  different  network  configurations  on 
estimates  of  the  errors  involved  in  locating  a radiating  source. 


? 

k 
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The  results  of  the  error  analysis  were  tested  by  using  statistical  distri- 
butions of  time-differences  (AtA  B)^»***»  etc.,  that  in  turn  yielded  a 
distribution  of  source  location  coordinates  and  signal  velocity  (<}><,,  0g,  v)^ 
For  given  means  and  variances  of  the  time  differences  their  normal  (Gaussian) 
distributions  were  obtained  by  taking  finite-sample  means  from  random  numbers 
for  which  routines  were  available  at  AFGL's  CDC  6600  digital  computer.  For 
certain  purposes,  rectangular  distributions  were  also  used. 
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In  order  to  illustrate  the  results  of  the  error  analysis  and  the  above 
tests,  coordinate -transformation  and  plotting  programs  were  developed, 
furthering  the  purpose  of  providing  a clear  visualization  of  network  performance 
and  testing  procedure.  The  initial  source  location  S (4>g,  0^)  was  placed  at 
the  center  of  a grid  system  whose  elements  were  made  to  represent  latitude  and 
longitude  error  boxes  expressed  by  Oq  , Oq  . These  displays,  which  comprise 
10  x io  standard  deviations  and  thus  exaggerate  the  smaller  error,  were  obtained 
by  rectilinear  projection  from  their  positions  on  the  curved  surface  of  the 
earth  onto  a plane  tangent  to  S,  such  that  great-circle  paths  through  S would 
appear  as  straight  lines. 

The  results  of  this  study  were  presented  to  Commission  VI  of  V.R.S.I. 
during  the  1973  International  IEEE/G-AP  Symposium  and  V.S.N.C./V.R.S. I.  Meeting 
in  Boulder,  Colorado. 
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Eigenvaiue-E-igtnvecto'i  PeXeAmination 


Initiator:  Mr.  G.  Borgiotti 

Problem  No:  4721  Project  No:  4600 


This  problem  involved  finding  the  eigenvalues  and  eigenvectors  of  the 
problem: 


xj  Vx)  /J  *4(0  « 


The  Initiator  wanted  accurate  values  for  the  largest  eigenvalues  and  their 
associated  eigenfunctions.  Towards  this  end  the  function  i|k(x)  was  represented 
by  a large  number  (N  up  to  2001)  of  equally  spaced  points  and  the  largest 
eigenvalues,  eigenfunctions  were  determined  as  follows: 


Set  up  an  equally  spaced  grid  of  X values  X. 


For  the  1st,  3rd,...  (2n+l)th  eigenfunctions  (j)  make  an  initial  guess 
of  a constant  for  the  values  of  ipj  (XjJ . For  the  2n d,  4th, ...  (2n)th  (j) 
eigenfunction  make  an  initial  guess  ^jCXi)=Xi 


3.  For  every  eigenfunction  except  the  first,  reset 


j-1 


'J'jtxp  = ^(X.)  - Z ^(xp  ("x—i ^ ) 

L-  X 


J[j  ^(5)  *l(5)  dC 

t— 


where  \|j  . (£)  4>T  d&  is  approximated  using  a trapazoidal  rule 
-i  j L 

4.  Store  ^(Xi)  in  H>S(X.) 

5.  Calculate  and  reset  (X^)  through 


i , sin  c(X.-C) 

VV  - ^ f-  l-r-T  ^ 


again  using  a trapazoidal  rule  approximation 
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6.  Repeat  step  3. 

7.  Using  the  trapazoidal  rule  approximation  calculate 


FKK1  = (X)  <PS(X)  dx 

and 

FK2  = ^ ^l  »jW)2  dx 

and 


FK21  = f_\  (^s(X))2  dx 


8.  Check  if  1 1 - p^Tpj^y  | < ERR  where  ERR  is  an  input  parameter  when 


this  condition  is  satisfied  X. 

J 


FKK1 

Fk2l~7T  ’I'j  (X)  is  finally  reset  to 


i 

r 


f- 


i 


l 

\ 

t 


♦j  W 


OnH) 

" FK2*4 


If  condition  8 is  satisfied  increment  j and  go  to  step  2. 

If  condition  8 is  not  satisfied  go  to  step  4. 

The  results  of  this  analysis  were  reasonably  accurate,  since  the  eigen- 
values checked  with  eigenvalues  found  by  an  independent  method.  The  program 
did  have  problems  when  c=8  since  as  c gets  larger  X2  and  X4  have  eigenvalues 
which  are  nearly  equal.  This  prevents  the  procedure  used  from  separating  the 
2°d  and  the  eigenfunction.  Consider  the  following: 

oo 

Let  F (X)  be  the  n^  guess  of  4/ . (X) . Then  F (X)  = E a \J/.(X)  where  if 
n J n £=1  e 3 

j is  odd  a2=a^=. . .=a2n=. . .=0  since  Fn(X)  is  even  and  if  j is  even  a^=aj= . . . = 

a2n+l=..*=  0 since  Fn(x)  is  odd.  Further  the  procedure  reset  F (X)  so  that 

al=a2=* • •aj_i=  °*  With  this  knowledge  the  function  Fn+i(Xj,  to  the  accuracy 

of  the  integration  satisfies  the  formula. 
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Fn*l<X>  * K ^ S Xe  VX> 

where  K is  a constant.  If  A^  » A^^  iterations  quickly  lead  to 

fN(X)  - Kx  *j(X) 

where  is  a constant  but  if  A^  ~ A^^  then  the  iteration  leads  to 

Fn(X)  - KjCa^W  . ajt2  *j.2W) 
where  K is  a constant. 

The  initial  guess  for  (x)  was  good  enough  so  that  a2  > 100  a^  so  the 
problem  was  not  as  serious  as  might  otherwise  have  been  the  case. 

This  program  has  been  run  with  c=4  and  c=8  for  the  first  5 and  8 eigen- 
values, eigenfunctions,  respectively. 


Coupled.  Mode  Propagation  oh  Hydro  magnetic.  Wave* 

Initiator:  Dr.  H.  Radoski 

Problem  No:  4729  Project  No:  7601 


There  are  two  fundamental  hydromagnetic  modes  called  the  poloidal  or 
isotropic  and  toroidal  or  guided  modes.  In  an  inhomogeneous  medium,  such  as 
the  earth's  magnetosphere,  these  modes  must  always  be  coupled.  Simple  geometric 
models  of  the  magnetosphere  have  been  developed  to  study  such  coupled  mode 
propagation.  The  case  to  be  treated  is  the  cylindrical  model.  No  analytic 
solutions  to  the  wave  equations  have  been  discovered.  It  is  felt  that  numerical 
solutions  should  afford  valuable  insight  into  the  development  in  time  and  space 
of  the  electric  and  magnetic  fields. 
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The  purpose  of  this  problem  is  to  perform  the  required  numerical  analysis 

and  computer  programming  necessary  to  calculate  and  plot  P(x,t),  T(x,t)  as 

functions  of  x between  x = 0 and  x * x for  a set  of  values  of  t between  t = 0 

o 

and  t * tMAx»  P(x,t),  T(x,t)  as  functions  of  t between  t * 0 and  t * t|^  for 
a set  of  values  of  x between  x * 0 and  x = xq,  such  that  the  following  partial 
differential  equations  are  satisfied: 


(32/3t2  + n2x'2  + m2)  T(x,t)  = m(3/3x)  P(x,t) 


[32/3t2  + n2x”2  - x"l(3/3x)  x(3/3x)]  P(x,t)  = -mx-1(3/3x)  xT(x,t) 


A sample  case  with  the  values  used  as  input  parameters  and  the  boundary 
conditions  that  were  satisfied  is  presented  below: 


m = 1,  n = 1,  xQ  = 1,  t^  = 20 


T(x,t  = 0)  = (3/3t)  P(x,t  = 0)  = (3/3t)  T(x,t  » 0)  = 0 
P(x,t  = o)  = JjCwx),  where  Jj(u))  * 0 

T(x  = o,t)  = P(x  =*  o.tj  = P(x  = Xo,t)  =*  o 

The  resulting  calculations  are  then  to  be  checked  by  the  following 
associated  conservation  theorem: 


$xdx  = constant 


4>  = G «■  I 


G = T(x,t)  + (BZ)2 
I = P(x,t)2  + (BR)2  + (BF)2 
(BZ)  = (m/x)  £ T(x,t)  dt  . 

(BR)  - (m/x)  P(x,t)  dt 
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(BF)  * [3/3x  P(x,t)  - mT(x,t)]dt 


EG  = jTX°  Gx  dx 
El  = J£X°  lx  dx 

Difference  methods8  were  applied  to  the  coupled  partial  differential 
equations  such  that  the  error  term  was  of  the  order  64,  where  6 represents  the 
difference  in  the  independent  variable  between  two  adjacent  grid  points.  This 
approach  was  unstable  in  that  the  calculated  solutions  became  divergent  at  the 
boundary  x*0.  Analysis  of  the  coupled  equations  and  intermediate  calculations 
seem  to  attribute  this  divergence  to  the  term  (1/x)  (3/3x)  P(x,t).  This 
approach  was  then  modified  by  re-evaluating  the  solutions  (P,T)  at  a particular 
time  (t)  for  distances  x=6,26, . . . ,ifi  from  a polynomial,  with  no  linear  term, 
fitted  to  the  above  solutions  at  distances  x=0,  (i+l)6,  (i+2)6,...  this 
modification  was  successful  in  preventing  the  divergence  and  yielding  possibly 
acceptable  results.  However,  this  later  procedure  was  discontinued  in  preference 
to  finding  a method  for  calculating  the  solutions  at  all  points  and  without 
such  smoothing,  nevertheless,  it  was  encouraging  to  find  such  simplification 
so  helpful. 

I 

It  was  next  decided  to  utilize  a predictor/ corrector  method  (1)  based  on 
the  following  correction  formula: 

P(x,t.)  = (P(x,t.+1)  + 6P(x,ti_1)  - 4P(x,t._2)  ♦ P(x,t._3))/4 

T(x,t.)  = (T(x,t.+1)  ♦ 6T(x,t._1)  - 4T(x,t._2)  ♦ T(x,t._3))/4 

This  formula  was  derived  by  equating  difference  approximations  for  the 
second  order  partial  with  respect  to  time  at  the  t^  ^ increment  in  time.  It 
has  an  associated  error  term  of  the  order  61*.  The  procedure  followed  was  to 
calculate  values  for  P and  T at  t^+j,  use  the  above  formulas  to  correct  the 
values  at  t^,  recalculate  the  values  at  t^+j,  compare  the  later  two  predicted 


1 Numerical  Methods  for  Scientists  and  Engineers,  R.W.  Hamming,  McGraw-Hill, 
1962,  p.  186-210. 
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values  at  t^+1>  and  should  the  comparison  be  adequate,  proceed  to  the  t^+2 
evaluation,  else  perform  the  iteration  of  correcting  the  values  at  t^,  etc. 

The  choice  of  correcting  the  ti  points  rather  than  the  t^+1  points  was  made 
in  view  of  having  two  successive  points  for  adjustment.  In  addition,  it  was 
preferred  to  provide  the  iteration  option  than  resort  to  a variable  time  step. 

The  validity  of  resulting  calculations  of  P and  T were  supported  by 
both  their  satisfying  the  conservation  theorem  and  their  repeatability  with  a 
change  in  the  grid  size  (6) . During  the  processing  for  various  combinations 
of  the  parameters  and  initial  conditions,  it  was  noted  that  as  the  parameter 
related  to  faster  damping  of  the  P wave  amplitude  approaches  10,  the  solutions 
diverged  at  the  boundary  x=l,  for  the  grid  size  6 = .001  in  both  space  and 
time  dimensions  as  time  (t)  approached  10. 


The  computer  core  and/or  time  requirements  seem  to  preclude  any  further 
significant  reduction  in  6 or  extended  use  of  the  iteration  method.  And  so, 
in  order  to  meaningfully  extend  the  range  of  processing,  it  seems  a variable 
grid  size  should  be  used  in  the  spacial  dimension,  especially  near  the  boundaries 
and  a re-analysis  of  these  coupled  partial  differential  equations  is  required 
in  order  to  find  a more  suitable  correction  formula. 


As  can  be  readily  appreciated,  a significant  computer  programming  effort 
was  required  to  develop  a structure  which  could  efficiently  store  and  retrieve 
data  for  both  calculation  and  presentation  purposes.  In  particular,  to  provide 
the  options  for  selecting  when  (time  tQ)  and  where  (distance  xq)  the  solutions 
P and  T are  to  be  outputted  in  plotted  form. 

VIZ. , letting  i (upper)  and  j be  indices  respectively  referring  to  time 
and  space  incrementing  in  the  difference  grid;  then,  excluding  points  near  a 
boundary  the  predicting  formula  used  for  P(x,t)  is  the  following: 


pr  ■ 2pj  - pi'1  - st2  pj  - (-pw  * spj.i  - 8pj.i  * 

(125x  Xj)  - (Pj+1  - 2Pj  + pj_j)/6x2 
♦ m£C-Tj+2  + 8Tj+i  - 81j_i  ♦ Tj_2)/126x  ♦ Tj/x.]  J 
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Pouka  SptcX/wJL  KnaJLyt>-U  Tide  Data 

Initiator:  Mr.  David  Anthony 

Problem  No:  4754  Project  No:  8607 

This  problem  involved  the  Power  Spectral  Analysis  (P.S.D.)  of  tide  data 
presented  as  a function  of  time.  The  data  was  initially  presented  to  us  on 
paper  tape  but  was  first  copied  onto  magnetic  tape  before  beginning  any 
analysis.  The  Maximum  Entropy  Method  was  used  to  determine  the  P.S.D.'s  and 
the  following  is  the  procedure  used: 

Initially,  the  data  was  read  from  the  magnetic  tape  and  edited.  This 
involved  determining  "bad"  or  missing  data  points  and  linearly  interpolting 
to  supply  "corrected"  points.  Occasionally,  the  editing  detected  large  gaps 
in  the  data,  in  which  case  the  data  was  processed  as  if  each  group  of  data 
came  from  a separate  paper  tape.  The  edited  data  was  then  placed  on  what  was 
called  the  "first  master"  data  tape.  It  should  be  noted  that  the  tide  data 
was  edited  and  stored  on  this  "first  master"  in  increasing  time  order. 

The  edited  data  was  then  re-read  and  an  interpolation  performed  between 
the  sets  of  data.  The  interpolation  was  performed  by  modifying  a theoretical 
tide  subroutine  supplied  by  the  problem  initiator.  Letting  T(t)  represent 
the  tide  value  associated  with  time  t,  then,  this  program  determines  values 
for  A,  B and  At  such  that  A+B*T(t+At)  is  a "best  fit"  in  a least  squares  sense 
of  the  data  on  either  side  of  the  gap  in  the  data.  The  program  then  finds 
values  of  A1  and  B1  such  that  A1+B1*T(t+At)  touches  the  data  point  immediately 
preceding  the  data  gap  and  touches  the  data  point  immediately  following  the 
data  gap.  The  values  of  the  function  A1+B1*T(tl+At)  at  the  times  at  which  data 
points  are  missing  are  the  interpolated  data  values. 

This  data  set  is  now  stored  on  the  "second  master"  tape,  from  which  a low- 
pass  numerical  filter1  is  applied  to  the  data.  The  pass-band  was  set  less  than 
* where,  Ny  is  the  original  Nyquist  Frequency  of  the  data  and  NDEC  is  a 
decimation  factor,  which  is  an  input  parameter  to  this  program. 


1 K.W.  Behannon  and  N.F.  Ness,  "The  Design  of  Numerical  Filters  for  Geomagnetic 
Data  Analysis",  NASA  Technical  Note,  NSDS-TND-33411. 
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The  data  is  then  decimated  by  NDEC,  and  M coefficients  (M4  is  pre- 
determined prior  to  execution)  determined  from  the  maximum  entropy  method 
(MEM)  are  applied  to  the  decimated  data.  (A  similar  approach  is  used  in 
the  solution  of  Problem  No.  4810,  described  later).  A Fast  Fourier  Transform 
is  now  used  to  determine  16385  P.S.D.  values  at  equally  spaced  frequencies, 
spaced  between  zero  and  the  new  Nyquist  Frequency. 


The  maximums  found  above  are  then  analyzed  by  finding  the  "exact"  frequency 

of  the  maximum  (i.e.,  the  peak  power  point).  This  is  accomplished  by  a search 

in  the  region  of  the  maximum  until  P„  , > .99  * P < P , where,  P is  the 

M-i  — m — m+1  m 

power  at  the  peak  power  point,  and  Pm  ^ (Pffi+1)  is  the  power  at  the  next  lower 
(higher)  frequency  studied.  In  addition,  the  area  (energy)  in  the  vicinity 
of  the  peak  is  calculated.  This  area  is  between  points  at  half  the  peak  power 
or  a minimum  of  a power  versus  frequency  curve,  whichever  definition  defines 
a smaller  integral. 

I 

The  results  showed  that  the  high  frequencies  agreed  with  theoretical  values 
to  within  .1%  but  the  low  frequencies  were  inaccurate.  The  energy  under  the 
peak  also  followed  the  expected  pattern.  We  expect  to  improve  the  accuracy  of 
results  with  an  improved  MEM. 

i 
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A LcueA  Beam  Slewing  Through  a.  TuAbuZent  Medium 

Initiator:  Dr.  R.  Fante 

Problem  No:  4756  Project  No:  2153 

When  a laser  beam  is  slewed  through  a turbulent  medium  the  spectrum  of 
the  intensity  and  phase  fluctuations  are  affected.  The  purpose  of  this  problem 
is  to  evaluate  the  integral  equations  which  describe  the  frequency  spectra  of 
the  log-amplitude  and  phase  fluctuations.  For  a beam  slewing  at  an  angular 
rate  tos  the  frequency  spectra  of  the  log-amplitude  and  phase  fluctuation  Qg 
is  known  to  be 


sin2  (ie=& 

2 1+a2 


cos2  (1^5  (i!2!i)) 


1+a 


where  ft 


radian  fluctuation  frequency, \/o  = ambient  atmospheric 


wind  speed,  L * path  length  in  turbulence,  k = signal  wavenumber,  y * u L/r  , 
2 so 

F *>  L/(kLg),  Lq  * outer  scale  size  of  the  turbulent  eddies. 


It  is  readily  noted  that  there  is  a singularity  at  the  lower  limit  of  the 
t integral.  That  is,  at 


t = ft1 2/(l+y£)2 

Analysis  of  the  integrand  reveals  it  is  of  the  form  1AJI.  Hence,  it  was  found 
acceptable  to  increment  this  limit  by  adding 


Lni2l/_«L_+p\ 

\l+Y£)2  / 


n/ie 

+ Fl  (RF) 


with  RF  representing  an  assumed  insignificant  fraction  of  the  integration. 
The  infinity  at  the  upper  limit  of  the  t integration  was  replaced  by 


This  was  determined  by  noting  that  the  exponential  term  of  the  integrand 
eventually  becomes  dominant. 

These  double  integrations  were  satisfactorily  evaluated  using  the  Adaptive 
Simpson  Rule! ' technique  developed  by  this  laboratory. 

1 A New  Adaptive  Simpson  Integration  Routine,  Neil  Grossbard,  Space  Data 

Analysis  Laboratory,  Boston  "College, "Chestnut  Kill,  Mass.  02167;  AFCRL-70-0504, 
Scientific  Report  No.  1,  September  1970. 


73 


A Functional  Determination  fio*  Fitting  Roentgen  Measurements 

Initiator:  Mr.  R.  Frederickson 

Problem  No:  4766  Project  No:  5621 


The  purpose  of  this  problem  is  to  determine  a functional  curve  fit  to 
data  which  represents  the  dependence  to  exposure  readings  in  Roentgens/ 
minute  (Ro)  on  the  radial  distance  from  the  collimator  center  line  to  the  ion 
probe  center. 

The  given  raw  data,  Ro,  is  first  corrected  by  the  factor  R/R1  which  is 
defined  as  follows: 


-1.25 

R = J ^ E N(E)  y(E)  dE 


-1.25 

R1  = J S(E)  E N(E)  y(E) 


N(E)  = nle" 


(,1<E<1.25) 


N(E)  = e (E=l . 25) 

(n  is  an  input  parameter) 

I = 1.25  [sin(1.37(E-.l))]l/3 

y (E)  = 4.505  x 10‘2  - 2.325  x io"2E  + 1.344*10_2E2 
-3. 728*10-3E3  - 3.64*10"3E_1  + 1.66*10"4E_2 


S(E)  = 8.733*10_1  + 4.999*10“2E  + 8.234*10"2E2 

-3. 178*10-2E3  + 3.053*10"2E-1  - 9.856*10"4E"2 
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Notice  that  N(E)  has  a "delta  function"  contribution  at  E=1.25.  Analysis  of 
the  above  integrals  reveal  that  it  is  satisfactory  to  utilize  Simpson's  Rule 
for  their  numerical  evaluation  using  the  continuous  form  for  N(E)  and  adding 


(i)  1.25  u(1.25)  exp(-N)  to  R. 

(ii)  1.25  S(1 . 25)  p(1.25)  exp(-N)  to  R1. 

The  procedure  chosen  to  determine  the  functional  curve  fit  to  the  data  is 
to  minimize  the  least  square  expression 


* = £ (MO  - F(t„,a,)) 


where,  M equals  the  number  of  data  points,  t is  the  independent  variable  (radial 
distance) , F the  mathematical  model  function,  and  a^  the  parameters  in  that 
function.  A rather  remarkable  curve-fit  to  this  data  is  obtained  upon  using 
the  following  mathematical  model 


F(t,a.)  = a + a exp(a  exp(a  ta1*)) 

J o 1 2 3 


A computer  program  incorporating  the  Space  Data  Analysis  Laboratory  version 
of  applying  the  Method  of  Levenberg  to  certain  nonlinear  problems  readily 
determines  appropriate  values  for  the  parameters  a^ . Included  in  the  output 
of  this  program  are  the  description  of  the  fitting  expression,  the  estimated 
values  for  the  parameters,  the  input  data,  the  fitted  evaluations,  the  cor- 
responding differences,  the  associated  variances,  the  standard  deviation  and 
plotting. 


Evaluation  and  Plotting  o£  Sett>  oi  double.  I ntegnali 

Initiator:  Dr.  G.  Borgiotti 

Problem  No:  4773  Project  No:  4600 


This  problem  involves  the  evaluation  and  plotting  of  sets  of  double 
integrals.  The  plots  are  generally  of  g(US  |sin9o)  versus  US  where 


g(USjsin0o)  = f(Xjsin0o)  e 


j(us1)x1 


, -jS  K sin20 

£(xJsin0o)  ■ f.  1 exP  UnT  V * — 7-1  C5J-2WX  5 ) 

a ’ y iii 


1 + -L  sin20 

I 2 a 

! y2 


E dp  + 1/2  exp(j(p+j)  nf^  - d 


given  j - V-T  , y=6.25,  sin0a=l/2,  Ka=314.16  and  sin  e=.08  and  the  values  of  S 
dp+1/2  are  input  parameters. 

In  order  to  calculate  f(x^|sin0o)  the  value  of  the  real  and  imaginary 
parts  of  f(x^|sin0o),  fr(xjsin0o),  fj(x  |sin0o)  are  calculated  separately. 
Thus,  let 


ARGi  = n£  v + 


S*  2 

1 + — sin  0a  and 

,,2 


Let  ARGd  = (p4)  tt(C  - 5i5|2.) 


Then 


i 


f (X  | sin0o)  = sin(ARGl)  Z dp+i  sin(ARG  ) 
1 p=-4  2 p 


+ cos(ARGl)  Z dp+-y  cos(ARG  ) 
p=-4  2 p 


and 


fj(X  | sin6o)  = cos(ARGl)  Z dp+y  sin (ARG  ) 
1 p=-4  p 


- sin(ARGl)  Z dp+y  cos(ARGp) 


p=-4 


The  integrals  for  fp(X  |sin6o)  and  f (X  |sin0o)  are  approximated  by  an 

K l I j 


adaptive  Simpson  Rule1  technique  developed  by  this  Laboratory. 


Values  of  f(X  |sin0o)  are  calculated  for  a given  value  of  sin0o  X = -1, 

1 1 1 
-1+AX;  -1+2AX,...  0,  AX,...  1 where  AX  = N set  through  an  input  parameter. 


These  values  of  f(X  |sin0o)  are  then  multiplied  by  weights,  as  if,  to  approximate 


-1 

Jx  f(Xjsin0o)  dXt 


by  Simpson's  Rule.  The  first  N+l  of  these  weighted  f(X  |sin0o)  (X  varies  from 
-1  to  0)  are  placed  in  the  first  N+l  locations  of  a complex  array  f*.  The  last 


N values  of  the  weighted  f(Xjsin0o)  (X  varies  from  AX  to  1)  are  placed  in 


location  4097-N  through  4096  of  f*.  Locations  N+2  through  4096-N  of  f*  are 
set  to  zero  and  a 4096  point  Fast  Fourier  Transform  (F.F.T.)  is  applied  to  f*. 
The  resultant  answers  contain  g(USjsin0O)  for  values  of  US  = where 


L varies  between  -2048  and  2048. 


A New  Adaptive  Simpson  Integration  Routine,  Neil  Grossbard,  Space  Data 
Analysis  Laboratory,  Boston  College,  Chestnut  Hill,  Mass.  02167;  AFCRL-70-0504, 
Scientific  Report  No.  1,  September  1970. 
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The  method  used  to  calculate  g(US  jsinOo)  is  a general  method  which 
can  be  used  to  get  sets  of  integrals  of  the  general  form 


X(s)  = / y(t)  ejts  dt  . 
a 


\ 


The  results  of  this  problem  were  accurate  plots  of  g(US  |sin0o).  One  plot 
for  each  of  a set  of  values  of  S and  0o. 

l 


A nalysis  ofa  MuZtt-Exponenttal  Luminescent  Decay  Rates 

Initiator:  Mr.  R.  O'Neil 

Problem  No:  4789  Project  No:  8658/CDNA 


When  optical  radiation  at  a given  wavelength  may  be  the  result  of  one 
or  more  specific  atomic  or  molecular  transitions,  the  magnitude  of  the  various 

contributing  emitters  may  be  established  by  analysis  of  the  time  dependent 

* 

optical  radiation.  The  analysis  is  simplified  if  the  number  of  contributing 
radiating  states  is  known  and/or  if  the  transition  probabilities  are  appreciably 
different.  Recovery  of  the  number  and  magnitude  of  the  luminescent  decays  as 
well  as  the  exponential  term  coefficients  will  then  be  further  analyzed  as  a 
function  of  molecular  collision  frequency  to  provide  reaction  rate  coefficients 
for  various  processes. 

The  purpose  of  this  problem  is  to  determine  the  number  and  magnitude  of 
luminescent  decay  rates  in  oxygen  and  oxygen-nitrogen  gas  mixtures  for  the 
case  where  the  decay  may  be  described  by  either  one,  two  or  three  exponential 
rates.  This  requires  developing  or  applying  existing  programs  to  various 
experimental  data  to  solve  an  expression  of  the  form: 

-x  t -x  t -x  t 

I = I e 1 + Ie  2 + Ie  3 +1 

2 3 4 
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1 


where  I is  given  for  various  values  of  t,  and  the  r,  and  their  standard 
deviations  are  to  be  determined. 

The  data  was  presented  in  two  forms:  (i)  tabulated,  (ii)  a data  tape 

representing  an  analog/digital  conversion.  The  first  form  could  be  readily 
processed;  whereas,  the  second,  which  was  output  by  the  AFGL  ac/dc  facility, 
required  both  unpacking  and  sorting.  The  unpacking  was  done  by  programs 
obtained  from  the  Boston  College  SDAL  Satellite  Group.  The  resulting  data 
tape  was  processed  such  that  the  data  representing  the  magnitude  of  the  lumines- 
cent decay  is  sorted  in  terms  of  time  identifying  the  beam  status  on/off  and 
then  placed  onto  another  tape. 

Data  sets  selected  from  this  later  tape  were  then  curve  fitted  by  the 
exponential  programs  developed  by  the  Boston  College  SDAL  Numerical  Analysis 
Section.  For  each  data  set,  the  curve  fitting  was  for  all  the  data,  some  of 
the  data  as  "averaged  data"  and  including  or  excluding  constraints  on  the  param- 
eter, in  particular,  the  parameter  representing  the  background  level. 

The  output  of  this  rather  substantial  processing  was  provided  in  both 
printed  and  plotted  form.  In  addition  to  plotting  on  a linear  or  regular  scale, 
an  option  was  also  included  for  displaying  the  results  using  logarithmic  scaling. 
This  was  accomplished  by  taking  the  logarithm  of  the  difference  between  the 
data  and  the  background  level  and  selecting  proper  intercepts,  such  that  the 
decay  rates  could  be  interpreted  as  the  slopes  of  straight  lines  which  correspond 
to  the  so-called  best  fittings.  A composite  plot  was  also  generated  which 
showed  the  dependency  of  the  slow  decay  rate  on  the  omission  of  data  pertaining 
to  a fast  decay  rate  as  a function  of  time. 

The  results  of  this  effort  were  presented  at  the  meeting  of  the  American 
Geophysical  Union  in  December,  1975,  at  San  Francisco,  California.  It  is  also 
intended  that  these  results  be  published  in  an  AFGL,  Scientific  Report. 
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Ele.eXAjon.-l on  MaXlnemaXieaZ  UodeXLLng  in  the.  lonatpkeAe. 


Initiator:  Dr.  John  Jasperse 

Problem  No:  4799  Project  No:  8627 


The  purpose  of  this  problem  is  to  find  a mathematical  model  to  explain  the 
behavior  of  electrons  and  ions  in  the  ionosphere.  Towards  this  end  many  models 
were  investigated. 

The  first  model  tried  assumes  the  ionosphere  contains  monoenergetic  photo- 
electrons. This  model  calculates  the  normalized  electron  velocity  distribution 
function  (Fnorm  (X))  and  the  electron  temperature  compared  to  the  photo- 


electron temperature  ( (T  /T  ) (y)).  The  equations  used  are  the  following: 

c pc 


F (X)  = 
norm 


i + (W  6 


{ 5 (X-X0)  + Z B 6 (X-(Xo-m)]  } 

m=l 


(T  /T  ) (y)  = + Z C (y)  [1-~mX°1] 

P 1 + y v^Xo”  m=lm  1 + y (X0-m)' 


+ cM(y)  [i-M  x0  1 ] 


.hereB  . S f [*•  ~ , 

m 1=1  »•  Vr/V1  + [Xo  - (£-1)]' 


C (y)  = " ( 

m n . \ - 1 

£-1  y + 


^i)jj j 

[X,  - (t-l)P 


Here  V = a n with  a the  recombination  coefficient  and 
r e 


ng  the  recombination  frequency. 


Vj  = /XT/Lx  with  Xi  the  threshhold  velocity  squared  and 
Li  the  mean  free  path. 

Xo  is  the  photo-electron  velocity  squared 
M is  an  arbitrary  input  parameter 


■ 


m 


In  addition  to  the  above, 
1 


I = 


M-l 
+ l B 


1 ♦ 


(Vj/Vy)  A7  m=l 


m 


1 ♦ (Vj/Vr)  (Xo-m) 


, ♦ B 

*5  m 


is  calculated.  I should  equal-  one,  if  the  equality  stands  up. 

The  second  model  tried  assumes  the  ionosphere  contains  exponentially 
distributed  energies  of  photo-electrons. 

The  solution  of  this  model  is  approached  as  n goes  to  infinity  of  the 
velocity  distribution  function. 


F (X) 

norm  v J 


( 


/xT  x0 


) 


’ sr_ 

} ♦ (vi/vr)  J 


sr  e (/f  -i) 


] 


(1  - exp  (-Xo"1))'1  exp  (-X/X0)  (X) 


(N) 


and  the  electron  temperature  compared  to  the  photo-electron  temperature 

(N)„  . vYT  rb 


T v"VT 
e pe 


Xo 


/ a*  F (N) 


norm 


Where  0 (X)  = 1 if  X >_  0;  * 0 otherwise 

B(N)  (X)  = [1  - exp  (-Xo-1)]  ♦ exp  (-NX0_1)  CN  (X) 

N-i 

+ [1  - exp  (-Xo1)]  E exp  (-nXo"1)  C (X) 

n=l  n 


Cn  (X)  = " , 
m=l 


it 

(x+m)  * 


Vr/V1  ♦ (X+m) 


the  input  data  is  the  same  as  the  first  model,  with  b the  upper  limit 

of  the  ratio  of  V./V  . 

l r 

The  integral  is  calculated  as  the  sum  of 

51  f\x  51  ,(«_  .51  fi%  51  ,00 


Xo 


norm 


Xo 


f 


norm 


with  each  integral  evaluated  using  a trapazoidal  rule. 
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I 


In  addition  the  following  quantities  are  calculated: 


,(N)  7 


£ ^“xT  _(N) 


Idx  ’ (X)  where  this  check  should 

/ 0 r-s-  norm  v 

o £.Y  A 


approach  one  as  N approaches  infinity 


F(o)  (X)  = f L_  ) 

norm  /x7  X0  1 + (Vj/V^  JT  0 (/)T-1) 


[1  - exp  (-Xo'1)] 


exp  (-X/Xq) 


I<°>  = 

(Z  * 

0 2 sr 

F(°) 

norm 

(X) 

T (0) 
e 

/XT  rb 

v/>C 

F(°) 

r 

pe 

a UA 

Xo  o 

~r 

norm 

and  the  electron  source  function 


(V  G ) (X)  = (2an  / (X0/x7)  ) /T  exp  (-X/X0) 

norm  e 

The  third  model  calculates  the  photo-electron  source  function 

So  (z)  = QpL1  n^z)  Ibl  exp  [A^z)]  where  nj  (z)  = naJ  exp  [(a-zj/h^ 

and  Aj  (z)  = hj  Qpal|sec  x|  [^(b)  - n^z)]. 

We  also  calculate  the  spherical  electron  flux 


Pe  ■(*) 


1 -?  f-s-l 

to  F.  (z)  I y(z) 

r ion  J L J 


[S0 (z)  - - S2  (z)  ] and 

2 


calculate  the  electron  density  n (z)  = (— )**  y(z)  p (z)  and  calculate  the 

e 4 e 


the  electron  temperature  Tg(z)  = 3.2991  (10 
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Here 


here 


a is  the  base  height  of  the  ionosphere 
b is  the  base  height  of  the  exosphere 


QpL1  is  the  photo-electron  ionization  cross-section 

Qpaj  is  the  photo- electron  absorption  cross-section 

nal  is  the  neutral  particle  density  at  a 

hj  is  the  scale  height 

is  the  photo-electron  flux  at  b 

X is  the  zenith  angle 

a is  the  recombination  coefficient 
r 

Fion(z)  is  the  value  of  the  fractional  ion  concentration  at  the  height 
S2  (z)  » 0 

y(z)  is  the  electron  velocity  at  the  height  z. 

y(z)  is  either  a set  of  input  parameters  or  if  the  temperature 

is  calculated  then: 


-i-  = JL  { 1 ♦ Z 1 C (X(z))  — .t1.'”.  Xo-^ , 

y(z)  y0  1 ♦ X(z)  SK  m*l  m 1 ♦ X(z)  (Xo-m)* 

♦ C M (X(z))  [1-MXo'1] 

where  X(z)  ■ V'  L'~\z)/ar  Fion(2)  ne(2)  and 


L‘‘(z)  - Q'  n'(a)  exp  ((a-zj/hp 


and  C (X(z)) 

in 


m 

n (• 
£-1 


[Xo  - (£-1)] 


X(z) -1  ♦ [Xo  - (£-1)] 

Yo  is  the  photo-electron  temperature 

Q'  is  the  inelastic  scattering  cross-section 

n'(a)  is  the  neutral  particle  density  at  a 

V'  is  the  threshhold  velocity 

Xo  is  the  square  of  the  photo-electron  velocity 

M is  an  input  parameter  (M  -*■  “) 
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The  fourth  model  solves  the  following  integral  equation 

pe(z)  = ri  + Kn  C1!!  + Fn  PeCz)  ♦ Fi2  CPeCz)*)) 

b 

where  Kn  = - J dz'  Ej  (—  | e'z/h  -e~z'/h|) 

2 a Le£ 

Ej  is  the  exponential  integral  of  the  first  kind 

Ij  = i(in)  exp  C — — (exp( - )-  exp  ( - )) 

[ VoLeJl  h n 

I11  = s<>  (z')  - ~ S,(z') 

2 £ 

. -V.2  y2(z') 

F11  = fLeZ  (z,)  + (S_1)  LiX  (z ')  e 1 (1  + V 

ira 

Fi2  = y2(z’)  Pion  (z') 

4 

further  (z»)  - IT*  exp  (-  ) 

h 

L:1  (Z-)  = L'1  exp  (-  -1'  ) 

h 

with  L’J  = Qe£  n(a)  exp  (2-  ) 

h 


= Q.  n(a)  exp(  — ) 


is  the  incident  flux  at  b 

Vo  is  the  cosine  of  the  angle  of  the  incident  flux 
f is  an  additive  fraction  in  the  equation 
h is  a scale  height 

n(a)  is  the  neutral  particle  density  at  a 


Qei  is  the  elastic  cross-section 


Vj  is  the  ionization  threshhold  velocity 
S is  the  mean  number  of  electrons  per  ionizing  collision 


To  solve  this  integral  equation  we  consider  a monotonically  increasing 
set  of  NVAL  zk  values  starting  at  a and  ending  at  b.  is  then  approximated, 
except  at  z = z ' , by 


. rALj-  (E,  -«!- 1 

K»1  2 L,t 


with  a,  - . and  = ^VALTI  ‘ b 


A (zk>  *7  (zk.l  - zk-l>  « z ' Z1 


the  term  in  the  sum  is  approximated,  using  a taylor  series  expansion,  by 


— E1  | e-z/h  -e'zk/h|)  A (zR)  as  z - Zj 


el 


1 Ej  (—  | e_z/h  -e"zk/h  | ) Af  (zk) 

2 Lefc 


♦ - Ej  (5l  | e~z/h  -e"zk/h|)  Ah  (zk) 
2 Lel 


where  Ej  (—  | e'z/h  -e'zk/h|)  At  (zk)  (i=f  or  h) 
Le* 


e.‘ 


= h (-E. , - e.  An  (5l  ) - I ( -i  )+  e'z/h  (2  - e x'<  . e x/‘) 


i/2  -ei/2 


L.i  2 2 


-e 


i/2 


. e cexP  (-2nz/hA  e.  - (1-e  i/2)  (1  + hi j 


n=l  2n  (2n)  J 


•)(  i 
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e.,,  , i/2 

♦ (1  - e l/2)  (1  ♦ 1 ~ e 


00  ,-(2n-rl)z/h  . 

- Z exP  ( ) 

n=2  (2n-l)  (2n-l)  ! 


-ei  /9  i i/2 

((1  - e 1/2)  (1  + i-® 


. ei/2 

+ (1  - e x/2)  ( 1 + ^ 


))  + 2 Z — (1  - e'Mei/2) 

2 n=l  nz 


here  Af  (zR)  = - (zR+1  - *k) 


(zk)  = “ (zl 


" 2v_i) 


Error  bounds  were  found  for  truncating  the  above  infinite  sums  and  the 
above  term  was  evaluated  to  a fractional  accuracy  of  .0001. 

With  this  approximation  for  NVAL  equations  were  formed,  one  for 
each  value  of  z^.  This  set  up  NVAL  equations  in  NVAL  unknowns.  We  solved 
this  equation  as  follows: 

After  setting  up  the  equation  NVAL  initial  guesses  for  the  values 
°f  Pe(zk)  are  made.  These  initial  guesses  are  either  input  data  or 


pe(V  = 


4 (So  (z.)  - S/2  S2  (z,)) 


a tt  F.  (z.) 
ion  v k' 


A "better"  guess 


is  then  calculated  as  described  below: 

For  each  equation  find  the  amount  (E(z))  for  which  the  equation  is 
not  satisfied.  Then  calculate  all  the  possible  derivatives 


3E(z) 


3p(zk) 


Then  form  NVAL  equations;  one  for  each  z value  considered. 


NVAL  3E  rzl 

z A z (z) 

k=l  dp  (z  ) 


-E  (z) 


”1 
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and  solve  the  NVAL  linear  equation  in  NVAL  Az^  unkowns  for  ' 


the  Az^  s.  Then  form  new  p C values 


PNEW  * pe^zk^  * a Az  calculate  the  new  ECz^J's 

^^NEW  (zk^  * Here  a *s  one  but  is  halved  until 

NVAL  2 NVAL 

E (E  (z  ))2  JVAL  f .2 

k=l  NfcW  K “b1  '^k'-'  . This  is  the  "beter"  guess  for 

p(z).  This  process  is  repeated  until,  judging  by  the  amount  of  error 
and  the  changes  in  p(z)  the  process  appears  to  have  converged.  As  in  the 
third  model  y(z)  is  either  an  input  parameter  or  is  calculated  as  indicated 
in  equation  (1).  In  addition  to  p (z)  and  y(z)  the  program  also  calculates 

k 

the  electron  concentration  n (z)  = (— ) y (z)  p (z),  and  the  electron  temperature 

4 e : ;.i 

12  * 

T (z)  = 3.2991  10“  / Cy  (z))  • If  y(z)  is  changed  the  program  does  each 

correction  of  Pe(z^)  after  changing  y(z). 

The  fifth  model  solves  the  following  integral  equations 


Pe  (z)  = Ij  * 

K11  «ll  + F11 

pe(z)r) 

+ K12  «21  + 

F21  pe(z)  J(z)  + 

F22  J(z)) 

and  J(z)  = I2 

+ K2i  (Xii  + Fn 

P(Z)  + i 

+ K22  (I21  + 

F21  p(2)  + F22  J 

(z)) 

where  Kj2  = — 

/ dz'  e2  c—  i 

A I 

e-2/h  -e-2 

2 a Le£ 

» i t < d,'  E,  |e-2/h  -e-z'/hl)  si*" 

2 * Le* 


(2) 


K22  .Ifiz-  E,  |e-2/h  -e-zVh  |) 


Je£ 
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where  sign  (X)  = 1 for  X > 0 and  -1  for  X < 0 and  where  I:  is 

— n 

the  exponential  integral  of  the  nth  kind 


Further  y0  exp 


fh 


do  L 


(exp  (-  - ) 


el. 


exp  (-  ^ )) 
h 


2a. 


21  = - — 3 y2(z')  F._  (z'); 


a on 


N 


F„  = -(1-f)  L-J  (z ' ) - E L'1  (z')f  (V.y(z') 


el  ^ J ~ £ “k  v“  3 M '■’V 


-l:1  (*')  f4  (v. y(z*)); 


where  f4  (x)  = — ±-  [x3  e“X*  + ix  e"x2  ♦ (i_erf  (X))] 

3/r  2 4 


l:1  (z*)  = l:1  , z'  . 

k k exp  ( ) 

h 


Lk  = \ n(a3  exP  W for  K = 1,  2,  N 


Here  N is  the  number  of  inelastic  cross-sections  used  in  the  analysis 
is  the  Kth  inelastic  cross-section 
is  the  Kth  threshhold  velocity 


K I 


This  set  of  integral  equations  is  solved  exactly  analogous  to  the 
fourth  method  but  2NVAL  equations  in  2NVAL  unknowns  are  formed.  The  additional 
unknowns  are  the  J(z^)'s  and  the  additional  equations  come  from  solving  (2). 

If  J(zk)  is  to  be  initialized  it  is  calculated  from  the  initialized  pe(z^) 
values  as  the  approximation  of 


J(zk)  - I2  + K21  C1!!  + F11  PeW  + Fl2(Pe(z3323  + K22  1 


21 
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J(zk)  we  calculate 


After  "solving"  for  pe(zk)  an^  the  particle  current  density 

H Jfz^ 


fz.  J = (-  j yfz.„)  ofzv)  the  electron  drift  velocity  U(zk)  * 
4 


W 


andTe(zk)  » 3.2991  10'1!/(y(zk)! 


From  the  mathematics  supplied  by  the  Initiator,  the  following  equation 
serves  as  a check: 


ttoi 


J(b)  - J(a)  = /bdz*  {So  (z’)  - (—  ) Fion  (z')  y2(z')  P2(z') 


-V?  y(z')2 

♦ (S-l)  L'1  (zr)  e (1  + V?  y2(z'))  p(z')  } 


This  equation  has  been  checked  with  the  integrals  calculated 
using  a trapazoidal  rule. 

Another  set  of  equations  derived  from  the  supplied  mathematics  should 


also  check: 


3 TTQ 

— J(z)  = So(z)  - (—  ) Fion  (z)  y2  (z)  p2  (z) 


3z 


(3) 


* (S-D  Li1  (z)  (1  + V?  y2  (z))  e 1 


-V2  y2(z) 


P(z) 


for  all  z values.  This  equation  is  checked  for  the  NVAL  z^ 
values  considered.  For  the  purpose 


— J(z. ) for  K t or  NVAL  is  approximated  by  fitting  a parabola  through 


3z 


the  points  zR  j,  zk  and  zk+J  and  then  finding  the  derivative  of  the  parabola 
at  the  point  zk.  The  derivative  at  z ^ and  z^^  is  approximated  as 


a J(Z-)  - J(z.) 

ST  ’ 


and 


z2  - 21 
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J 


3 ^ J(ZNVAL5  ' J(ZNVAL-1) 

7"  JUNVALJ  ' 

oZ  Z. 


"NVAL  ' ZNVAL-1 

The  program  may  reset  y(z)  as  in  the  third  model  or  by  finding 
y(z)  values  which  satisfy  (3).  When  equation  (3)  is  used,  equation  (3) 
is  not  used  for  checking.  In  any  case  the  program  now  solves  for  the  electron 
temperature  Te(zk)  = 3.2991*  10“12/(y  (z^))2. 


The  sixth  model  solves  the  following  integral  equations  in  a 
manner  exactly  comparable  to  that  used  in  the  fourth  and  fifth  models. 

pe(z)  = - — J2(z)/pe(z)  + Ij  + Kn  Ijj  + K12  I21  + K13  I31 
31T 


* K11  F11  pe(z*> 

+ K11  F12  pe 

(z«) 

♦ K12  F22  J(z') 

+ <K11F13  + 

K13 

k12  21  ^e' 

T2  frn 


* f,4  * K13  F32> 


J(2)  = l2  * K2x  xn  + K22  *21  * K23  *3X 


* *21  F11  »e<z'>  * K21  F12  * K22  F21  pe‘z'>  J<z'> 

* k22  f22  j(z')  . (k21  f13  . k23  f31)  Jz(2-) 

* <K21  F14  * "23  F32> 


Further  I21  = Sj  (z1)  = 0 

I3I  a S2  ( z * ) = 0 
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F,J  - ar  /*(*')  Fi0„  (••)  FiOT  («•) 

A N • X f 7 M 

F..  « (— ) { (4-3f)  L~J  (z')  ♦ Z a L.  (z')  e vk  y 1 

14  Sir  61  k-1  K 

-V*  y2  (z*) 

[2  ♦ 2V2  y2  (z*)  ♦ Vj  y1*  (z')]  * L71  (z‘)  e 

[S  ♦ 3 ♦ (S  ♦ 3)  V2  y2  (z’)  ♦ 2S  V*  yk  (z*)]  } 

2a 


*31-  -~T  y2  (Z,)  Fion  <*'>* 

N -V2  v2  i z ' 1 

F32  " {1-^  Le£  (z'>  * E Lk*  (z,)  e k 

32  15ir  * k-1  K 

-V2  y2  (z') 

(1  ♦ Vj  y2  (z')  ♦ j V*  y4  (z1))  ♦ Ljl  (z')  e 

(1  ♦ V2  y2  Cz')  ♦ - y-  (z*))  } 

l 2 i 


After  "solving"  for  p#  (zfc)  and  the  particle  current  density 
J(zfc)  we  calculate 


r H 

Ne  (*k)  ■ F ) Y Uk)  P (*k); 

U (zk5  " J (zk)/Ne  (zkJ  ; 

Te  (zk)  - 3.2991  • 10'12/  (y  (zk))2 

and  the  spherical  electron  flux 

pe  (zk}  " p (zk}  { 1 + (“  ) j2  tzk)/p2  (zk5 
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The  integral  and  derivative  checks  explained  in  model  five 
also  are  used  in  model  six  as  may  be  the  change  in  y(zk)  from  equations  (1) 
or  (3) . 

When  y(zk)  is  incremented  from  (3)  then  the  derivative  check  is 
not  used.  Further  when  y(zk)  is  increment  then  Te(zk),  Ng(zk)  and  U(zk) 

are  recalculated  after  the  changes  before  a new  iteration  on  p (z, ) and 

0 1C 

J(zk)  is  begun. 

The  final  model  tried  is  the  same  as  model  six  except  pe(zk)  is 
not  calculated  and  the  change  in  y(zk)  which  occurs  is  the  implementation  of 
the  following  equation.  This  equation  is  approximated  using  methods  previously 
discussed. 


(-a. 

15ir 

) ifi2!  = i 

P(2) 

+ K31 

N11 

♦ 

K32 

N21 

+ 

K33 

N31 

+ K31 

1H11  ^ <*') 

N 

+ E 
k=l 

[»!! 

v; 

/ 

(z‘) 

+ 

Hk 

H12 

V2  y2'  (z*) 

♦ *4si 

-V?  y2  (z«) 

i 

, -1 

y2  (z*) 

e * 

+ Vll 

vi 

(z') 

+ 

H12  Vi 

y2 

(z’) 

* H13>e 

1 

+ K32 

{ H2j  y2  (z*) 

~N 

+ E 
k=l 

[»21 

Vk  : 

v3  ( 

z’) 

♦ 

k 

"22 

Vk  y (z’))e 

-Vj.y2(z') 

* H23 

(1  - erf  (Vk  y 

’ ((z*))]  + 

<H21 

V?  y 

3 (z 

+ h22  vi  y 

(z')) 

exp  (- 

■V2  y2'  (z’))  + 

H23 

(1  - 

erf 

<vi 

y (: 

Z’)))  > 

* K33 

{ m31  y2  (z* ) 

N 

+ E 
k=l 

t»31 

K 

y4 

(z*) 

♦ 

Hk 

"32 

V2  y2  (z’) 

*H33 

-v"2  y2(z') 
] e k 

+ I»31Vi 

y4 

Cz’) 

♦ 

H1 

"32 

V2 
! vi 

y2  (z*)  ♦ 

H33^ 

-V?  y2  (z') 
e 1 } 


A 

b 

■ 

y • 
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4 


■i 


4(S+3) 


Lj;1  (z*)  J2  (z')/p  (z') 

) l:1  (z»)  j2  (z')/p  (z') 

) L?1  ,z')  J2  (z')/p  (z') 

) l:1  (z«)  j2  (z')/p  vz') 

The  results  of  the  methods  discussed  varied  from  case  to  case.  The 
two  methods,  model  S and  model  6,  judged  to  be  the  closest  approximations 
to  the  actual  ionosphere  were  close  approximations  of  each  other.  For  this 
reason  the  investigator  judged  that  the  particular  approach  which  was 
programmed,  has  been  pushed  as  far  as  reasonable  at  the  present  time.  The 
final  result,  model  6,  did  not  conform  to  the  experimental  results  in  the 
ionosphere.  The  investigator  thinks  this  is  because  the  magnetic  field 
was  not  incorporated  in  his  model.  The  investigator  intends  to  add 
the  magnetic  field  and  return  for  programing  services  for  the  improved 
model. 
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1.  BACKGROUND 


This  request  for  numerical  analysis  and  scientific  computer  programming 
support  involves  the  derivation  and  application  of  certain  product,  reciprocal 
and  natural  logarithm  algorithms  to  the  exact  evaluations  of  derivatives  of 
given  functions. 

The  primary  function,  (see  mathematics,  J function),  is  defined 

and  used  by  Jasperse*  in  the  calculation  of  Atomic  Coulomb  Integrals  (ACI) , 
*V"on"on'ono”  T^ese  integrals  arise  in  studying  the  physical  properties  of 
ozone,  an  important  constituent  of  the  earth's  atmosphere. 

In  the  original  analysis/programming  support  request,  the  J function  was 
specified  by  the  Initiator  as  a function  of  4 variables,  (a,  3,  Y,  6),  where 
the  computational  requirement  was  to  calculate  "high  order"  derivatives  of  a 
function  of  four  variables.  Certain  derivatives  of  the  J function  subsequently 
would  be  applied  to  computations  of  L integrals.  However,  preliminary  inspec- 
tion of  the  J function  showed  that  it  was  necessary  to  treat  four  dummy 
constants,  (a,b,c,d),  as  variables  when  taking  partial  derivatives.  Therefore, 
the  problem  was  immediately  expanded  to  . requirement  for  evaluations  of  high 
order  derivatives  of  a function  of  eight  variables,  a considerably  more  com- 
plex undertaking. 

i 

It  should  be  noted  that  the  J function  contains  logarithmic  terms,  an 
additional  problem  to  be  dealt  with. 

It  should  further  be  noted  that  previous  attempts  by  the  Initiator  to 
compute  ACI  (to  which  the  J function  is  applied)  by  direct  integration  techniques 
had  met  with  limited  success;  low  order  computations  had  been  generated  on  a 
computer,  but  the  higher  order  calculations  were  beyond  the  scope  of  the 
techniques  and  equipment  available  at  that  time.  These  low  order  calculations 
were  presented  to  be  used  in  a study  of  accuracy  at  the  time  when  the  evaluated 
J function  would  be  applied  to  ACI  calculations. 


2.  PROCEDURE 

The  synoptic  approach  taken  to  the  problem  is  outlined  below.  Some  of 
the  steps  included  were  not  known  at  the  onset  of  the  effort,  but  were  made 
necessary  by  difficulties  encountered  as  the  analytical  investigation  proceeded. 
Obviously,  most  of  the  effort  outlined  below  progressed  in  a parallel  fashion 
rather  than  in  specifically  sequential  steps. 

Evaluations  of  J Function  Derivatives  and  Application  to  Calculations 
of  L Integrals. 

1.  Define  a basic  method  of  evaluation  of  J function  derivatives  by 
applying  product  and  reciprocal  forms  of  algorithms. 

2 

2.  Expand  existing  algorithms  (Product,  reciprocal)  from  functions  of 
1 variable  to  functions  of  8 variables. 

3.  Prove  the  validity  of  the  above. 

4.  Determine  suitable  approximations  for  terms  such  as 

W ^%oS'a~r  2(f^rATelj-J  Which  can  reasonably  be  evaluated  on 
the  computer  by  the  methods  product  and  reciprocal  algorithms. 

5.  Apply  evaluated  J function  to  computations  of  certain  low  order 
ACI.  It  is  necessary  to  have  prepared  compatible  programs  or  sub- 
routines to  evaluate  Coulomb  integrals,  Sturmian  functions  and 
Gegenbauer  Polynomials.  Generate  study  of  accuracy  of  results 
compared  to  previous  calculations  by  direct  integration  technique. 

6.  Determine  an  efficient  computer  index  storage  scheme  to  allow 
evaluation  of  large  numbers  of  partial  derivatives  of  high  order 
J functions. 

7.  Derive  algorithm  for  exact  evaluation  of  derivatives  of  terms  such 

,,,  .-Polynomial  1(8  Variables), 
as*  LN  ‘•Polynomial  2 fit'  Variables) } ' 

8.  Prove  validity  of  the  above. 

9.  Apply  logarithmic  derivative  algorithm  to  appropriate  portion  of  J 
function;  compare  accuracy  of  results  to  approxiamtion  techniques. 
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10.  Apply  logarithmic  derivative  algorithm  to  low  order  ACI  calculations; 


compare  accuracy  of  results  to  previous  study. 

11.  Compute  high  order  ACI. 

3.  MATHEMATICS 

3.1  The  following  function  (J  Function)  is  given  by  the  Initiator1 
[p.  17,  Eq.  19]: 

yySciE  _ rr^ 

2(fga2b2c262-ega2b2d2Y2-fb2c2d2a2+ea2c2d2B2) 

X Jin 

where  a,  6,  Y,  <5  are  expansion  parameters;  e,f  are  particle  mass  ratios; 

a,  b,  c,  d are  dummy  variables  introduced  in  a derivative  substitution 
technique. 

3.2  The  Atomic  Coulomb  Integral  (L  integral)  is  given1 [p.  17,  Eq.  20]: 


(fabciS+abdY+acdB)  (gabc6»bcda+eacd$) 


dy+abc6+bcda) (gabdy+acc 


iY  5 o B _ jjY  ijO  x/d^  8 ci  B jY6d8)a  k c j 

V'on'Wono  V»o  n"o  Vo  no  A|V"n"n*n  J /a’D,Q'° 


where  N are  normalization  factors  of  the  Sturmian  function. 

I 

3.3  The  differential  operator,  D,  is  defined1  [p.  17,  Eq.  21]: 


I 


_ i n-P  P 

n_A  BP  a •»  a •* 

x Z s;  - (-2-)  (-2-) 

P *o  no  3B2  3b2 

4 


3.4  Derivative  Algorithms  for  Numerical  Evaluation 

The  method  taken  to  evaluate  partial  derivatives  of  the  J function  with 
respect  to  (a,B,Y»^,a,b,c,d)  consists  of  being  able  to  define  evaluations  of 
partial  derivatives  (to  as  high  an  order  as  required  by  the  full  evaluation) 
of  each  component  term  of  the  overall  expression,  and  then  being  able  to 
algebraically  combine  the  partial  derivatives  of  those  terms  in  a manner 
which  forms  the  evaluation  of  the  partial  derivatives  of  the  sum,  product, 
reciprocal  or  natural  logarithm  of  the  given  component  terms.  By  defining  all 
partial  derivatives  of  all  component  terms  of  the  J function,  and  by  using 
a high  speed  computer  to  evaluate  the  expressions  and  manipulate  the  combining 
of  the  expressions,  the  evaluations  of  the  partial  derivatives  of  the  full  J 
function  can  eventually  be  obtained. 

Obviously,  algorithms  are  required,  based  principally  upon  Leibniz's 
Rule,  to  form  the  evaluations  of  partial  derivatives  of  Products  of  terms  of 
polynomial  functions  of  eight  variables,  reciprocal  of  such  a term  and  the 
natural  logarithm  of  such  a term. 


The  algorithms  used  to  perform  these  operations  are  given  below. 


3.4,1  Derivatives  of  a Product  of  Functions  of  More  Than  1 Variable 


I , I ... 

v i 2 


D'  Yz  (f(x)g(x))  = Z Z 


H V2  ••7IlVI2V"  iX,i2  •** 


X ,x  .., 
1 2 


(f(x)). 


i =0  1,-0  °X  »X2  *•• 


I -i  ,1  -i  ... 

V 1 2 2 (g(x)) 

A • A • • • 

1 2 


where  D^f  (f(x))  represents  the  derivative  of  f(x)  with  respect  to  variable  X.. 

A 1 


r ft  T-  i ‘ifr'A  t.ii  ' - .... 


I 


3.4.2  Derivatives  of  a Reciprocal  of  a Function  of  More  than  1 Variable 


i ,i  ...  V1  V1  •••  /i \ /e  \ ••• 

- (.t).. 


i -l,l  -l  ... 

(1/fW)  D 1 2 (f(x))/f(x) 

A * A • • • 

1 2 


3.4.3  Derivatives  of  the  Natural  Logarithm  of  a Function  of  More  Than  1 


variable 


I ,1  ,1  ...  I*. 

,12  3 N 


I 1 1 I - 

i 2 3 ...  n 


v1  * x3  " j£n(f(x))|  = E E E E I1!)#1*!/1* 

VV  s *•  N i =o  i =0  i =0  •••  i =o 


1 2 3 


IM-l\  i ,i  ,i  . ..  iM 

N V.xV  •••  x >fMl 

L 1M  / 12*3  N 


I -i  I -i  I -i  ...  I„-iM 

v.x\; «■  N N 

12  3 N 


3.5  The  expression  selected  to  approximate  £n(x)  terms  in  the  initial 

3 

evaluations  of  the  J function  consisted  of: 


*nx  = 2 igr  + i®r>  + •••]  ^ 


Derivatives  of  this  expression  were  initially  obtained  using  the  product 
and  reciprocal  algorithms  defined  above. 
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3.6  An  example  of  introducing  substitute  derivative  expressions  with  y 
variables  to  evaluate  an  integral  expression  is  given  below. 


Given 


0*-^ 


S“  x2e~x  dx  = 


Introduce  dummy  variable  a 


Take  of  each  side 


Rewrite 


Take  ^ 
3a 


3 r°°  .-ax2 j 3 1 fir 

is  ■i  e dx  -sriVs 


•6  s|  e-axJdx  = - i 

3a  4 3/2 

a 


r x2e-ax  dx 


4 V/2 


Evaluating  at  a=l 


r"xV<‘>x2 . ‘^”.1  = 

0 4 Li  4 
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3.7  Derivative  Variable  Transformations 


Derivatives  were  evaluated  with  respect  to  first  order  variables  but  were 
required  in  the  differential  operator,  D,  with  respect  to  second  order  variables. 
Therefore,  the  following  relations  were  applied  to  each  term  of  the  Differential 
operator: 


i r 3 i1  - 1 3 

ax2  3X  ax1 

2.  (J-)2  = J_  (Jl)  . I -L 

ax2  4X2  ax2  ax1 


3.8  First  Order  Derivative  Expressions 


Expressions  for  certain  first-order  derivatives  of  the  J function  were 
derived  and  are  presented  below: 

Given: 


1 

fg<52  - egy2  + eg2  - fa2 


iln 


f ff6+v+g)  (gd+ot+eg) . 
l(ey+6+oO  (gy+S+fa) 


Rewrite : 


2__jY6ag  _ (fg52_egy2+eg2-fa2) 
ft2 


[4n(f6+y+g)  + S.n(g6+a+eg)  - 
jln(ey+6+a)  + £n(gy+g+foO  ] 


3.8.1 


J.  W.  [(-1) (fgd2-egy2+eg2-fa2)  * (2fg6)][in(  ) ♦ £n(  ) - *n(  ) - to(  )] 


♦ (fg<52-egy  +eg  -fot  ) tfd+y+g  + gd+a+eg  (ey+d+a) ^ 


3.8.2 


^ ^ [(-l)(-2)(2fg6)(fg62-egY2+eB2-fcx2)  * (-2egY)][iln  ♦ to  - £n  - in] 

* [(-l)(2fg6)(fg62-egY2+e62-fa2)"2][7^  - —J-  - -Jfogfi 

* [(-D(fg62-egY2+e02-fa2)  * (-2e*y)J “ 5yW 

* [(fg62-egY2*eB2-fa2)  *][ — LLf HK  ] 

(fS+y+B)2  (ey+S+a)2 


3.8.3 


2 3 3 3 3J  _ r(24»16)e2f2g2aeY<S1ro_  (f6+y*B)  (gS+q+eB), 

^2  3a  38  3?  37  " ^ D5  (eY+'fi+Qi)  (gY+B+fo)  3 

+ r48e2fg2gY<S1  r 1 1 f i 

1 d4  j lgo+a+e$  ey+6+a  gY+3+faJ 

f48ef2g2ot6Yi  r 1 . e 1 ! 

*■  q4  ■*  ‘■fS+y+B  g '6+a+eB  ~ gy+B+fa 


iiS£^-nlsi 


+y+B  ” eY+(5+a  " gY+p+ 


f48e2fga$Y1 > f g _ 1 , 

1 d4  j Lf6+Y+8  g6+a+eg  eY+°+aJ 


+ [3efg2y6j  j- e_ 


(g6+a+e0)2  (gY+3+fa): 

♦ [8e.f&e.Y][ £ ♦ £& .] 

D3  (ey+6+a)2  (gY+3+fct)2 

+ [8£.2jgy]  [ L 


(gS+a+eB)2  (ey+S+a) 
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. [ £ i_ 

D3  (gy+B+fa)2  (f6n*B)2 


. [Mgl]  t ®£ + I ] 

D3  (gfi+a+eB)2  (f6+y+B)2 


j8efaB_j  e 


-] 


(ey+6+a)2  (f6+y+8)s 

♦ [2£gi][ U& i + [1®J£L]  r 2eg , 

D2  (gy+B+fa)3  D2  (g6+a+eB)3 

♦ [2|e][ 2® _]  ♦ [¥±)[ if j 

D2  (ey+6+a) 3 D2  (f«+y+6)3 


where  D = fg6z-egy2+eB2-fa2 


4.  VERIFICATION  OF  SUBROUTINES  PRODUCT,  RECIPROCAL 

Suitable  polynomial  functions  were  selected  which  allowed  evaluation  of 
derivatives  of  their  products  and  reciprocals  by  analytical  methods.  As  it  is 
unrealistic  to  analytically  evaluate  derivatives  of  products  or  reciprocals 
of  expressions  of  functions  of  eight  variables,  functions  of  three  variables 
were  used.  This  was  felt  to  be  sufficient  to  show  that  some  method  for  dealing 
with  more  than  one  variable  is  either  valid  or  not  valid.  The  algorithm  derived 
should  later  be  able  to  be  expanded  by  symmetry  to  include  any  number  of  subse- 
quent variables.  In  addition,  the  particular  computer  software,  CDC  6600 
FTN  3.4.3,  was  found  to  be  bounded  in  subscripting  capability  to  three  indices. 

As  a minimum  of  time  was  desired  to  be  spent  on  the  programming  of  supportive 
test  functions  for  numerical  evaluations,  this  factor  was  taken  into  account 
when  initially  deciding  to  limit  the  testing  effort  to  functions  of  three 
variables. 

Computational  results  showed  that  the  derived  product  and  reciprocal 
algorithms  for  functions  of  three  variables  were  valid  when  compared  to  analytici 
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evaluations  involving  the  same  processes.  As  an  additional  test,  use  was  made 
of  the  property  that  "*■  0 where  K is  any  constant.  Considering  that  any 
indexing  scheme  for  control  of  storage  of  large  multidimensional  arrays  can 
readily  become  complex,  it  was  felt  that  a proof  based  upon  this  property 
would  lend  considerable  confidence  to  the  validity  of  the  computer  programming 
of  the  algorithms  into  general  storage  subroutine  form.  As  a guide  to  the 
Scientific  Programmer,  it  should  be  noted  that,  with  respect  to  all  of  the 
derived  algorithms,  the  value  of  the  differential  is  stored  in  an  array  whose 
indices  can  lead  to  the  order  of  differentiation. 

Therefore,  the  following  test  was  devised: 

Given:  f = f(x,y,z),  g = g(x,y,z)  and  all  derivatives  of  f and  g through 

3n 

3xn3yn3zn 

1.  Evaluate  all  f and  g at  some  (x  ,y  ,z  ),  and  store  in  multidimensional 

0 0 o 

arrays  F,G 

2.  Form  array  P = F*G  using  product  algorithm  subroutine 

3.  Form  array  R = using  reciprocal  algorithm  subroutine 

4.  Form  array  Z = R»P  using  product  algorithm  subroutine 

The  contents  of  array  Z will  then  be 

Z(l,l.l)  = 1. 

all  other  Z(i,j,k)  = 0. 

This  test  not  only  indicates  validity  of  the  algorithm  and  computer  sub- 
routines, but  also  indicates  any  level  of  roundoff  which  can  be  attributed  to 
the  methods  (in  conjunction  with  the  specific  functions  and  data  values  selected) 
for  the  steps  outlined. 


5.  PROGRAM 


' Ll,  .wu  - 


Functional  Description 

The  logical  computational  procedure  for  the  action  taken  by  the  computer 
program  can  be  broken  into  a series  of  sequential  steps.  Consider  that  the 
function  to  be  evaluated, J,  can  be  simplified  to  the  following  expression: 


A = Q(i)  An(Y)  - Q(^)  *n(j££) 


where  A,X,Y,M,N,0,P  are  polynomial  functions  of  eight  variables, 
and  Q is  a real  constant. 

Then  the  program  logical  procedure  for  numerical  evaluations  of  the  J 
function  and  application  to  AC1  can  consist  of  the  following: 


dn(0),  d"(P) 


dn(0*P) 


.n,  1 


1.  Calculate  and  load  into  arrays 

2.  Calculate  derivatives  of  product 

3.  Calculate  derivatives  of  reciprocal  : d“(Q~p) 

4.  Calculate  and  load  into  array 

5.  Calculate  derivatives  of  product 

6.  Calculate  and  load  into  array 

7.  Calculate  derivatives  of  product  : dn ^ (*0]*N 

8.  Calculate  and  load  into  array 

9.  Calculate  derivatives  of  recirprocal  : d“(i-) 

10.  Define  constant  Q 


: dn(M) 


AC^)  CM)] 


: dn(N) 


d"(X) 


jn,l. 
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11.  Calculate  derivatives  of  reciprocal  ; dn(^-) 

12.  Calculate  derivatives  of  In  term  ; dn[4n(Y)] 


1 


13.  Multiply  derivative  array  by  constant  Q 

14.  Calculate  derivatives  of  product  ; dn  [(Q)  £n(Y) • (~)  ] = dn(A) 

15.  Select  derivative  terms  to  apply  to  AC I calculations,  and  generate 
family  of  integrals. 

In  addition  to  several  test  sequences  of  certain  variables  which  demonstrated 
characteristic  properties  of  the  J function,  two  specific  data  sets  were  applied 
to  the  finished  program,  and  are  given  below  as  a matter  of  reference. 


VAR 

SET1 

SET2 

a 

2.31 

2.31 

8 

2.31 

1.14 

Y 

2.54555 

2.54555 

6 

.7093 

2.54555 

a 

1. 

1. 

b 

1. 

1. 

c 

1. 

1. 

d 

1. 

1. 

e 

.000136 

.999864 

f 

.000136 

.5 

Computational  Results  printed  by  the  program  include  a list  of  the  input  data 
parameters,  a list  of  all  partial  derivatives  of  the  numerically  evaluated 
differential  operator,  and  a list  of  values  of  the  associated  family  of  Atomic 
Coulomb  Integrals. 

6.  VERIFICATION  OF  J DERIVATIVE  EVALUATIONS 


To  test  the  evaluations  of  the  J function,  certain  low  order  terms  were 
determined  analytically  and  hand  evaluated  at  selective  data  points.  These 
included: 
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5a  5FJ 

,3  .2  3 Tooa$ 

l3oJ  3F  J 


3 ,3  Tooa0 

3aW  J 

,3  .2  r3  -w2  Tooa0 

(S3T>  (3F}  J 


evaluated  at  the  points  a*b*c=d=l 

e»f*l 

y«6»0 

a»1.0 

0-1.01 


In  addition,  the  analytical  expressions  for  the  first  order  derivatives  of  the 
J function  up  to  J^01^  were  determined  and  coded  as  a supplementary 

computer  program.  This  was  extremely  useful  in  debugging  the  initial  stages 
of  the  principal  program  to  evaluate  the  general  J expression,  in  that  computa- 
tions were  compared  for  validity  and  degree  of  accuracy.  Also,  insight  as  to 
the  behavior  of  the  component  terms  of  the  derivatives  of  the  J function  was 
gained  at  this  time,  which  was  to  prove  valuable  at  a later  time,  during  the 
study  of  computational  accuracy  of  higher  order  ACI. 

Other  analytical  studies  focused  upon  the  probable  cancellation  of  terms 
in  the  denominator  of  the  J function,  under  certain  conditions  of  equality  of 
pairs  of  terms,  as  suggested  by  the  initiator  on  the  basis  of  known  physical 
relationships  of  the  parameters.  These  include  (at  a=b=c=d=l)  a=0  and  Y=6. 

While  certain  terms  dropped  out,  no  obviously  significant  pattern  of  cancela- 
tion was  observed.  This  study  became  more  significant  during  the  study  of  the 
accuracy  of  the  ACI  evaluations. 

7.  COMPUTER  CORE  STORAGE 

In  programming  the  initial  studies  of  the  J function  it  became  obvious  that  ol 
the  twovstandard  principal  computer  difficulties  to  be  overcome  (central 
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processing  time/ core  storage)  one,  C.P.  time,  would  be  of  considerably  less 
impedance.  This  was  due  primarily  to  the  ease  and  speed  of  computation  that 
the  product  and  (most  particularly)  the  reciprocal  algorithms  allowed.  Even 
when  extrapolated  from  functions  of  3 variables  to  functions  of  8 variables, 
it  was  obvious  that  although  still  large  compared  to  other  major  efforts 
typically  submitted  to  the  machine,  the  C.P.  time  would  still  be  within  some 
realistically  manageable  range. 

The  core  storage  problem,  however,  when  dealing  with  the  derivatives  of 
a function  of  eight  variables,  immediately  attracts  one's  attention.  The 
following  chart  may  help  to  explain  why.  Consider  that  one  has  only  a function 
of  two  variables.  Then,  letting  F = f(x,y)  the  storage  required  to  hold  all 
of  the  evaluations  of  partial  derivatives  to  a given  order,  N,  can  be  graphically 
presented  as: 


Figure  1 - Partial  Derivative  Storage  Format 


Notice  that  the  location  containing  any  

3 ym  Sx11 


holds  the  value  of  the 


derivative  evaluated  at  some  given  (x,y).  One  sees  readily  that  the  space 
required  for  a complete  order  of  partial  derivatives  can  be  determined  by: 

Storage  = (order  of  derivative  + l)(No‘  of  variables) 

A storage  requirement  chart  can  quickly  show  the  limitations  imposed  upon  an 
array  size  in  terms  of  number  of  variables  and  order  of  derivative. 


NO.  \ 

OF  \ORDER 
VAR.  \ 


4 

16 

64 

256 

1024 

4096 

16284 

65536 


5 

25 

125 

625 

3125 

15625 

78125 

390625 


5 

6 

6 

7 

36 

49 

216 

343 

1296 

2401 

7776 

16807 

46656 

117649 

279936 

823543 

1679616 

S 764801 

Figure  2 - Core  Storage  Requirements  - 

No.  of  Variables  Vs.  Order  of  Derivative 


When  dealing  with  a function  of  eight  derivatives,  one  sees  immediately 

xd 

that  direct  core  is  no  longer  available  certainly  by  the  3 order  derivative, 
even  on  larger  machines.  One  should  also  bear  in  mind  that  the  storage 
requirements  given  above  refer  only  to  one  array.  Analysis  of  various  means 
of  evaluating  the  given  J function  showed  that  a minimum  of  4 such  arrays  would 
be  necessary.  This  meant  that  the  analysis  of  L integrals  would  be  limited  to 
approximately  (n* ' ,=n"  =n' =n=2)  if  the  entire  J function  would  be  calculated 
with  the  contents  kept  in  core. 


Otherwise,  several  options  are  open  to  the  programmer/ analyst: 


1.  Investigate  and  apply  external  core  storage  capability,  such  as 
magnetic  tape  or  disk.  (Virtual  Memory) 

2.  Break  analysis  into  smaller  steps  keeping  only  results  required  by 
further  analysis. 

3.  Determine  those  specific  parts  of  the  J function  required  for  evalua- 
ting the  L integral  (not  all  partial  derivative  components  of  the 

J function  are  used  by  the  differential  operator,  D),  and  find  a 
means  to  limit  the  calculations  to  those  levels. 

4.  Continue  to  solve  the  J function  analytically,  such  that  when 
eventually  applying  the  product  and  reciprocal  subroutines  the  numeri- 
cal analysis  starts  at  a higher  level. 

5.  Transform  the  J function  to  allow  derivatives  with  respect  to  a , 8 , 
Y2,  62,  rather  than  a,  8,  Y,  6. 

6.  Solve  the  J function  (or  the  L integral,  since  that  is  the  function 
to  which  the  J function  is  eventually  to  be  applied)  analytically, 
either  in  terms  of  some  variable  or,  preferably,  variables,  thereby 
significantly  reducing  the  overall  number  of  variables  to  be  kept 
in  storage  when  evaluating  derivative  expressions. 


8.  ATOMIC  COULOMB  INTEGRAL  CALCULATIONS 

As  shown  in  the  mathematics  section  of  this  report,  an  Atomic  Coulomb 
Integral  (L  integral)  can  be  calculated  by  applying  the  differential  operator, 
D,  to  the  derivatives  of  the  J function,  in  conjunction  with  the  appropriate 
normalization  factors  of  the  Sturmian  functions,  N. 

Those  AC I which  could  be  calculated  without  requiring  extensive  external 
computer  memory  subroutines  were  programmed  to  allow  evaluation  of  the  overall 
procedure  and  to  test  those  routines  required  for  storage  of  derivatives  of 
functions  of  eight  variables. 

Under  these  conditions,  all  ACI  from  order  L to  order  L were 
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successfully  calculated.  Comparison  of  some  of  the  lower  order  integrals  to 
results  previously  calculated  t other  techniques  indicated  that  the  accuracy 
of  the  newer  calculations  was  improved.  To  show  this  effect  more  clearly, 
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and  to  indicate  the  overall  degree  of  accuracy  of  the  derivative  technique , 
an  additional  version  of  the  program  was  generated  using  double  precision 
computation,  thereby  allowing  the  accuracy  of  each  term  of  the  L integral  to 
be  defined  when  compared  to  the  single  precision  computations. 

In  this  study  of  accuracy,  a pattern  of  decreasing  precision  was  observed 
which  eventually  limits  the  evaluation  of  the  L ingegral.  To  ensure  that 
such  loss  of  precision  was  not  a result  of  programming  errors,  a roundoff 
error  analysis  was  applied  to  both  the  single  and  double  precision  versions  of 
the  program.  The  minimum  detectable  delta  was  applied  to  each  input  parameter 
to  allow  observation  of  the  effect  of  a small  change  in  input  with  respect 
to  computational  output  over  the  full  range  of  evaluations.  The  results  of 
this  test  supported  the  conclusions  that  the  overall  precision  of  higher  order 
AC I computation  is  eventually  limited.  Whether  this  is  an  inherent  charac- 
teristic of  the  expression  being  evaluated  or  a function  of  the  applied  derivative 
technique  was  not  determine  at  this  time. 

' 

9.  CONCLUSIONS 
————— 

In  reviewing  the  results  of  the  application  of  derivative  algorithms  to 
a high  speed  computer  to  allow  evaluations  of  partial  derivatives  of  the  J 
function  and  the  numerical  computations  of  certain  ACI,  it  is  seen  that 
eventually  other  techniques  are  required  to  evaluate  higher  order  ACI.  As 
noted  in  the  study  of  computer  core  storage  requirements,  the  applications  of 
several  other  methods  of  combinations  of  numerical  analysis/scientific  pro- 
gramming techniques  could  be  investigated.  In  addition,  the  specific  reason 
for  the  eventual  loss  of  precision  in  this  analysis  should  be  determined. 

Of  the  options  available  from  the  above  methods,  that  most  likely  to  yield 
accurate  results  in  the  evaluations  of  significantly  higher  ACI  would  appear 
to  be  of  some  form  involving  item  6;  that  is,  an  analysis  which  reduces  the 
number  of  variables  involved  in  the  evaluation  of  the  J function  or  the  L 
integral,  yet  which  retains  a format  whereby  the  power  of  the  derivative 
algorithms  derived  herein  can  still  be  realistically  applied  to  the  computer 
to  allow  numerical  solutions  to  the  problem. 
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Under  this  consideration,  a method  outlined  by  Calabi  has  been  briefly 
reviewed,  where  an  outline  is  formulated  which  shows  that  ACI  can  be  solved 
in  closed  form  according  to  evaluations  of  sums  of  derivatives  of  functions 
from  which  selective  poles  are  removed,  as  typically  resulting  from  Cauchy's 
Residue  Theorem  as  applied  to  closed  loop  integration.  It  would  be  necessary 
to  investigate  whether  such  residues  could  be  evaluated  by  use  of  the  deriva- 
tive algorithms  already  derived  and  applied  to  the  J function  of  this  problem. 

In  addition,  it  would  probably  be  necessary  to  derive  a new  algorithm,  com- 
patible with  those  already  applied,  which  would  allow  evaluation  of  derivatives 
of  certain  functions  of  the  form: 

^-3-  |f(x,y)[  y = <Kx)]  x = x 

3x  3yJ  0 

One  advantage  of  being  able  to  apply  these  derivative  techniques  to  the 
Residue  Theorem  lies  in  the  fact  that  the  number  of  variables  over  which 
differentiation  takes  place  can  be  reduced  from  eight  to  possibly  two,  thereby 
1 exponentially  reducing  the  amount  of  core  required  and  immediately  allowing 
evaluations  of  much  higher  ACI. 

An  additional  advantage  would  appear  to  be  that  the  ACI  analysis  can  be 
carried  out  in  its  general  form,  that  is,  not  subjected  to  certain  constraints 
(subscripts  i,ll  = 0)  imposed  by  application  of  the  J function. 

On  the  surface,  it  appears  that  these  results  could  be  accomplished, 
thereby  allowing  the  initiator  to  calculate  full  families  of  ACI  which  could 
then  be  applied  to  the  solutions  of  higher  order  problems,  such  as  the  application 
of  a matrix  of  ACI  in  the  determination  of  the  bound  energy  states  of  certain 
molecules. 
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fou/UeA  AnafrpiA  oh  Wind  Senior  Data 

Initiator:  Dr.  E.  Dewan 

Problem  No:  4806  Project  No:  6687 

This  problem  involved  Fourier  Analyses  of  wind  sensor  data.  The  process 
used  in  this  analysis  is  to  find  the  average  Power  Spectral  Densities  (P.S.D.) 
of  many  consecutive  (8192  points)  sets  of  input  data.  The  averaging  procedure 
should  improve  the  stability  of  the  results. 

In  order  to  improve  the  results  of  the  analysis,  the  data  is  pre- 
whitened before  finding  its  P.S.D.  and  then  post-darkened  after  the  P.S.D. 
is  calculated.  Pre-whitening  is  useful  in  order  to  minimize  leakage  of  Fourier 
power  from  the  very  large  low  frequency  power  to  the  smaller  high  frequency 
power.  The  pre-whitening  used  is  to  analyze  first  differences  of  the  data.  The 

effect  of  using  the  first  differences  is  to  multiply  the  original  P.S.D.  by 
Fit 

2-2  cos  where  F is  the  frequency  in  cycles  per  time  unit  and  Ny  is  the 

nyquist  frequency.  Post -darkening  is  accomplished  by  dividing  the  resultant 
Ftt 

P.S.D.  by  2-2  cos  (— ) to  remove  the  effects  of  pre-whitening. 

The  results  of  this  problem  is  a P.S.D.  which  decreases  rapidly  with 
increasing  frequency  until  about  2/3  the  nyquist  frequency  and  then  begins  to 
level  off  with  frequency.  These  results  were  not  unexplainable  and  have  been 
used  by  the  investigator. 
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FouxLeA  Analytic  ofi  VoppleA.  Data 


Initiator:  Dr.  Kurt  Toman 

Problem  No:  4810  Project  No:  5631 

This  problem  involved  Fourier  analyzing  Doppler  and  various  control 
data.  The  Initiator  wished  to  study  how  the  Power  Spectral  Density  (P.S.D.) 
of  his  data  varied  with  time.  Towards  this  end  short  segments  of  the  data 
were  analyzed  for  their  (P.S.D.) . The  data  segments  analyzed  were  for  a 
particular  time  duration  (T) . Generally,  all  data  segments  for  which  the 
time  duration  (T)  of  the  data  was  available  were  analyzed  and  the  resultant 
P.S.D.  versus  freauency  was  plotted  for  each  possible  starting  time.  This 
set  of  plots  will  be  referred  to  as  the  isometric  plots. 

Originally,  linear  methods  of  Fourier  analysis  were  used  to  analyze 
the  data;  but  in  this  most  recent  problem,  a method  known  as  the  "maximum 
entropy  method"  (M.E.M.)  has  been  used.  M.E.M.  was  originally  suggested 
by  Burg  (1967) . In  the  analysis  for  this  problem  the  method  used  is  that 
described  by  N.  Andersen1 2.  Recently,  it  has  been  suggested  that  the 
algorithm  used  by  Andersen  is  inaccurate  (Dr.  Paul  F.  Fougere  of  AFCRL 
private  communication)  and  a "new"  algorithm  is  to  be  tried.  However,  under  our 
existing  contract  (F19628-73-C-0136)  the  Andersen  algorithm  has  already 
been  explored.  The  Andersen  algorithm  was  independently  derived  and  pro- 
grammed by  Dr.  Raj an  Varad^  of  this  laboratory  and  his  subroutine  was 
modified  and  used  in  the  studies.  The  flow  chart  derived  by  Andersen 
follows: 


1 On  the  Calculation  of  Filter  Coefficients  for  Maximum  Spectral  Analyses  by 
N.  Andersen  in  Geophysics  Vol.  39,  No.  1 (February  1972)  p.  69-72. 

2 

Data  Processing  With  Different  Techniques  for  Cross-Power  Spectra  by 
W.  Pfister,  G.S.  Sales  and  R.  Varad.  Environmental  Research  Papers  506  AFCRL 
Report  PR  75-0194. 
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t - l.M-l 

aa(t)  = 

a(t) 

t « l.N-M 

bl  (t)  » 

bl (t)  - aa(m-l)*b2(t) 

. b2(t)  - 

b2(t+l)  - aafm-l)*blft*n 

p * o 
t - 1,N 
P ■ p+X(t)**2 
P(o)  - P/N 
■ ■ 1 

bl(l)  - Xd) 
b2(N-l)  - X(N) 
t - 2,N-1 
bl (t)  * b2(t-l) 


om  * den  * 

t = 1 ,N-M 

nom  - 

nom 

+ bl (t)*b2(t) 

den  = 

den 

♦ bl (t)**2  + b2(t)**2 

a(m)  = 

■ 2*nom/den 

p(m)  = 

■ p(m-l)*(l-a(m)**2) 

t*l,m-l 

a(t) 


aa(t)  - a(m)*aa(m-t) 


From  a set  of  N equally  spaced  values  (x(t)),  this  algorithm  finds  M 
coefficients  (aj,)  and  the  square  residuals  of  the  linear  filter  (P^)  such 
that  for  a time  spacing  At  the  P.S.D.  for  any  frequency  P(f)  between  0 and 
the  nyquist  frequency  is. 


P(f)  = 


i1  - Z^V5 

m=l 


-2iri*f*m*At  1 2 


In  the  program,  written  for  this  problem,  P(f)  is  initially  calculated 
for  a grid  of  f values  (f^).  This  grid  is  usually  made  up  of  equally  spaced 

f values.  The  equal  spacing  is  modified,  however,  to  assure  plotting  at 

| 

frequencies  specified  by  the  Initiator.  The  results  are  modified  whenever 
a maximum  is  found  in  the  P(fj)  values  and  the  maximum  P(fj)  occurs  where 
3 < j < N-2.  In  this  case,  the  values  of  P(f._2),  P(f._j),  P(f..),  P(*i+1) 
P(fj+2)  are  modified  by  using  the  average  values  found  by  applying  a trape- 
zoidal rule  to  the  P(f)  between  f^  2 1 f 1 ^j+2*  '^lus » let 
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nr^r . pot ; tj2  < f < f(i '-1 


*1  * fi  1 i * f.r 


Pfr^T-  POT ; -i-~2-;  i 


f . , ♦ f j f . ♦ f . . , 

Ftfp-  - POT;  -i^—i  < f < -2_j-itl 


fi  ♦ f 


f,  . ♦ f. 


PPJ^T-  POT  ; 1 < f < ^ 


fTT^T-POT;  fj>2 


and  then  Pffj)  is  reset  to 


P<fj-2> 


P(fj_2)  (fj.2  ~ ^j-3^  * P(fj-2J  (fj-l  " fj-2j 


V^-3 


p<Vi} 


Pffj) 


mp 


P<W 


P(f  ) . P(fj+2J  (fj*3  ~ ^2j  * P^fj*2^  (fj»2  ~ fJ+l} 
j*2  *j+3  *j+l 


where  f0  - fj  and  fN+1  - f„ 

This  correction  is  necessary  because  the  maxima  entropy  Method  Makes 
the  P.S.D.  very  peaked,  such  that  the  amplitude  of  the  peak  is  not  a "good" 
Measure  and  that  is  what  the  correction  attempts  to  simulate.  This  idea  was 
gleaned  fron  reference  (3)  below. 


3 A Conparison  of  Power  Spectral  Estinates  and  Applications  of  the  Maxi mum 
Entropy  Method  by  Henry  R.  Radoski,  Paul  F.  Fougere  and  Edward  J.  Zawalick 
in  Journal  of  Geophysics  Research  Vol.  80,  No.  4,  February  1,  197S,  p.  619-625. 


119 


The  actual  program , which  will  be  submitted  to  the  SUYA  computer 
program  library  upon  completion  of  this  problem  proceeds  as  follows: 

a.  Linearly  interpolate  the  input  data  to  assure  equally  spaced  data. 

i | 

b.  Take  a first  difference  of  the  input  data.  The  program  processes  • ' 

these  first  differences.  This  procedure  has  been  found, 
empirically,  to  filter  out  a red-noise  behavior  of  the  initial 

data.  | 

i ' 

c.  Find  a set  of  coefficients  from  the  M.E.M.  for  the  entire  data 

I j 

set.  Calculate  the  P.S.D.  and  bar-graph  the  logarithm  of  the 
P.S.D.  at  a set  of  frequencies  g^.  The  g^'s  are  usually  equally 
spaced  with  a spacing  equal  to  a fraction  (input  data)  of  the 
resolution  frequency  found  for  a linear  Fourier  Analysis  of  the 
data.  The  spacing  of  the  g^'s  are  varied  to  have  the  plotting 
frequencies  specified  by  the  Initiator. 

d.  Calculate  and  plot  the  "isometric"  plots  P(f^)  described  above. 

e . Repeat  c . 

f.  Calculate  and  plot  the  results  of  a Kolmogoroff -Smirnov  Test  on 
the  results  of  c.  This  shows  whether  the  data  could  be  explained 

as  due  to  random  noise.  j 

! \ 

From  the  results  of  parts  c.  and  f.  above,  it  has  been  assumed  that 
the  procedure  using  the  entire  data  set  gives  very  stable,  and  probably 
accurate  P.S.D.'s.  The  results  of  d.  are  more  suspect  but  do  seem  to  allow 
identification  of  the  presence  of  power  at  particular  frequencies  and  how 
the  P.S.D.'s  behave  as  a function  of  time. 

To  test  the  results  of  the  program,  several  test  runs  were  made.  These 
runs  were  made  on  data  consisting  of  one  or  more  frequencies  plus  some 
random  noise.  The  random  noise  was  necessary  because  the  M.E.M.  breaks  down 
when  the  noise  level  is  zero.  The  test  results  show  that  the  M.E.M.  is 


better  than  a linear  Fourier  Analysis  in  predicting  the  frequency  of  an 
input  signal.  The  M.E.M.,  however,  did  show  soae  error  in  the  frequency 
found  and  some  double  peaks  where  single  peaks  should  exist.  It  is  hoped 
that  the  "correction"  suggested  by  Dr.  Paul  Fougere  will  provide  sore 
accurate  answers.  Towards  this  end,  a prograa  is  being  written  to  test 
Dr.  Fougere* s "correction". 


Analy&iA  Vapoun.  TxcuZ  Data.  F nom  PhotogKapiu 

Initiator:  Gordon  Best 

Problem  No:  4818  Project  No:  7635 


This  problem  involved  the  analysis  of  photographic  data  representing  a 
vapour  trail  released  by  a rocket.  The  task  was  to  determine  how  the  vapour 
trail  behaved  with  time  and  thereby  infer  the  velocity  versus  altitude  profile. 
The  difficulty  arose  from  the  fact  that  the  altitude  of  the  portions  of  the 
vapour  trail  seen  in  the  photographs  were  unknown. 

Attempts  were  made  to  fit  the  data  to  velocity  versus  altitude  models 
with  unknown  parameters.  These  attempts  were  unsuccessful  as  the  results 
indicated  that  a very  "good"  guess  to  the  solution  was  necessary  before  the 
parameters  could  be  determined.  A more  successful  procedure  eventually  was 
instituted. 

The  method  tried,  under  an  earlier  problem,  involved  finding  a wind 
velocity  versus  altitude  curve  for  each  azimuth  and  elevation  reading  determined 
from  the  photographs.  These  wind  velocity  versus  altitude  curves  are  the  locus 
of  points  which  would  lead  to  the  particular  azimuth  and  elevation  reading. 

Thus  let  Vj  j(A)  be  the  wind  velocity  of  the  Ith  photograph,  taken  at 
tne  Jfh  site  (assuming  altitude  A)  [There  are  three  sites  (J  values)];  then, 
for  a particular  azimuth  and  elevation  reading  say  M(A)  the  program  searches 
for  the  minimum  of  min  |VL  M(A)  - Vj  j(A) | where  M^J  and  with  this  restriction, 
we  search  over  all  possible  values  of  I,  J and  A.  The  program  indicates  the 
values  of  I,  J,  A and  ^(A)  determined  by  this  procedure.  The  point  VL  M(A),  A 
may  be  one  point  of  the  desired  velocity  versus  altitude  profile.  The  results 
of  this  analysis  gave  a rough  indication  of  the  actual  wind  velocity  as  a 
function  of  altitude. 


The  procedure  here  was  refined  as  follows: 

To  find  a minimum  of  min  C CVL  m(a)  " vi  j(A))2  * (VL  " Vs  t^A^2 
♦ (vr  r(A)  - V (A))2) . Where  MfJ,  Mj<t  and  J^t,  thus  we  consider  one  point 

1 f J S ^ u 

from  each  site.  Again,  with  these  restrictions  we  search  over  all  possible 
values  of  I,  J,  s,  t and  A.  After  finding  I,  J,  s,  t and  A by  this  procedure 
the  value  of 

VL.M<A>  * VI,jW  + Vs,t(A> 

y— 


is  considered  to  be  a possible  value  of  the  wind  velocity  at  altitude  A.  This 
procedure  leads  to  a more  refined  estimate  of  the  wind  velocity  as  a function 
of  altitude. 

At  the  present  time,  the  preceding  method  is  being  refined  to  interpolate 
between  azimuth  and  elevation  points.  We  hope  to  further  refine  the  wind 
velocity  profile. 

In  order  to  implement  the  following  procedure,  we  need  to  find  vl.m<a>- 
Thus,  let  Xr  be  the  position  of  the  rocket  at  time  tR.  Also  let  AZ,  EL,  t be 
the  azimuth, elevation  and  time  of  the  reading  L,M.  Further,  let  A be  the 
altitude  of  the  rocket  at  time  tR.  Then  if  L,M  is  due  to  a vapour  trail  at 
altitude  A,  its  coordinates  (in  rectangular  coordinates)  are  calculated  as 
follows: 

Let  $ be  the  geodetic  latitude  of  the  site 
A be  the  longitude  of  the  site. 

Then  from  (page  J,1  a rectangular  coordinate  system  of  this  site  can  be 
expressed  as 


1 Application  of  Vector  and  Matrix  Methods  to  Triangulation  of  Chemical 

Releases  in  the  Upper  Atmosphere  by  Antonio  F.  Quesada,  AFCRL  Environmental 
Research  Paper  No.:  351. 
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P = X 


[•--»  .k] 

LVl  - t2  sin2$  -1 


cos<|)  cosA 


,;[•  a • hi 
L\1  - e2  sin2$  J 


cos$  sinA 


, C T a(l-e2)  , "| 

♦ z ■ — ♦ h si 

. yjl  - e2  sin$ 


sin<J> 


where  if  the  shape  of  the  earth  is  assumed  to  be  an  ellipse 
a is  the  semi -major  axis  (equatorial  radius) 


b is  the  polar  radius 
e is  the  eccentricity 
h is  the  altitude  of  the  site 


and  therefore 


a|  * - — • a ~ * h|  sin<fr cosX  ♦ 

LVl-e^sin2^  J 


ae 


2 ~ r 1 5 sin * cos  * 
(l-e2sin2«j»)1,:> 


cosX 


dx 

dl 


f>/l-e2sin2<fj 


cos<p  sinA 


= sin<|»  cosX 


% * “ 1 ~ 8 - + hi  sin<j>  sinA  + 

(_Vl-e2sin24>  j 


3L£*  2 

cos  <j>  sin<|>  sinA 


(l-e2sin2$) 


ft- 


fefer] 


cos<f>  cosA 


^ = cos<()  sinA 
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JO  JO  is  Jr* 


• r ■ y.LV.LL.  ♦ hi  cos#  ♦ **  .CVf  J sin*#  cos# 

J l-c*sin# 

k w«n  of  AZ,BL  defines  • vector  in  the  direction  (PA  B) 

V . 1 

»A>£  • k (»U(EL» 

4 (cos(AZ)  cos(BL)) 

♦X  (-sin(AZ)  cos(BL))  £ 

t' 

A A 

#,  and  X ere  vector*  of  unit  length. 


i trensferuntlons  of  h,  #,  X into  X,  Y,  Z (vector*  of  unit  length) 
e*  is 

e sin(BL)  § ♦ (cos(AZ)  cos  (EL))  § ♦ (-sin(AZ)  cos(BL)) 

■ sin(BL)  J ♦ (cos (AZ)  cos(BL))  J ♦ (-sU(AZ)  cos  (EL))  jf 

Z,|  - sin(BL)  ^ ♦ (cos(AZ)  cos(BL))  ^ ♦ (-sin(AZ)  cos(BL))  § 


whore,  in  general,  J (a»x  or  jr  or  »,  b«h  or  # or  X) 

/■' 


12S 


mm 


a * ,da/  fda  2 

J?  * V * C3F 


AAA 


Therefore,  we  can  calculate  the  position  of  the  site  (Xg,  Yg,  Zg)  and 
the  possible  points  at  which  an  azimuth,  elevation  reading  point 


X * X ♦ BX 
p s r 

A A A 

Y «=  Y ♦ BY 

p s r 

A Ak  /V 

Zn  « Zc  + BZ„ 
p s r 


where  B is  a parameter. 

Next  we  consider  what  value  of  B would  explain  an  altitude  A for  the 
vapour  trail.  If  the  vapour  trail  is  released  at  altitude  A,  it  is  assumed 
to  be  blown  by  constant  winds  (for  the  time  interval  considered)  such  that  it 
stays  at  the  same  altitude  above  the  earth.  Thus,  the  position  of  the  vapour 
trail  is  on  sin  el  ipse  which  can  be  described  as 


A A ^ 

r = Z b sin  s + XC  + yD 
r 2 


C2  + D2  = a2  cos2 (s  ) 
r '2 


where  ar  = a+A2,  br  = b+A  and  s is  arbitrary.  Then  we  must  have  from  (1)  and 

(2) 

b.  sin(s  ) = Z + BZ„  ( 


aj  cos2 (s  ) « (Ys  + BYR)2  ♦ (Zs  + BZr) 


(4) 


Assume  we  guess  s *s*  B*B* 

2 


Then 


E = b sin(s*)  - Z - B*Z, 
x r ' 2 s F 
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— i - b cos(s*) 
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* " _ZR 


3E 


— §■  • -2a2  sin(s*)  cos(s*) 
3s  r 2 2 


3E 


3B 


I?  - -2  XR(Xs  ♦ “ Xjj)  - 2 Yr(Ys  ♦ B'YR) 


and  using  Newton-Raphson's  method  we  can  converge  to  a solution  for  B and  s^. 

£ 

Me  also  know  the  original  position  of  the  rocket  so  s * tan-1 
b_  1 

( — tan  A ) where  A is  the  latitude  of  the  rocket  at  altitude  A. 
ftj,  s s 

Then. assume  constant  latitudinal  winds  K (A)  at  altitude  A or 

x 


Cx(A)  = 
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1 ♦ — — - cos*  (s)  ds 


and  since  b ~ ar 
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1 r 1 

(S, 

- sin(o^)  cosCo^jl 


s ♦ sinfo  ) cosfs  ) 
1 2 2 


Then  since  we  assume  the  longitudinal  winds  are  constant  K (A)  at  altitude  A 

2 


^ b sin  s = K (A) 
ds  r 2V  * 


or  from 


dX  _ dX  dt 
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(A  -A  ) 
v p r 


a,  K.IA1  s2 
b k‘(A) 

1 c 


. b2-;.2  . b2-a2 

1,  r 1,  r r 


log  (.tan  (-j))  (1  ♦ + jl— — ~)  cos  (s) 


where  Ap,  Ar  is  the  longitude  of  the  vapour  trail  at  time  t,  rocket  at  time  tR 
with 


. „ - 1 ,Ys  + B\ 

Ap  ' tan  CXo  + BXr) 


Let 


S = (log(tan(^.)) 


b!-a_2 


log(tan(4)))  (1  + \ ^-L) 
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, b2-a2 
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♦ (COS(S  ) - 

^ _ 2 2 


cos(s^)) 


Then  finally 


b (A  -A  ) K (A) 
r P r i 


\W  ' a.  s 


Using  this  analysis  many  velocity  versus  altitude  curves  were  calculated, 
plotted  and  compared.  The  results  for  one  group  of  data  lead  to  an  estimate 
of  a wind  velocity  versus  altitude  profile.  The  results  of  another  set  of 
data  seem  to  indicate  some  inconsistency  in  the  data. 

This  work  is  continuing. 
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